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ABSTRACT 


Chaos metrics are examined as a tool to analyze 
atmospheric three-dimensional dispersion models at the 
individual particle rather than the aggregate level. These 
include the self-affine fractal dimension, D,, Shannon 
entropy, S, and Lyapunov exponent, A. Intercomparison or these 
metrics is first performed with the one-dimensional logistics 
difference and the two-dimensional Henon systems of equations. 
The fractal dimension and Shannon entropy are then measured as 
a function of the inverse Monin-Obukhov length (1/L) for two 
three-dimensional Lagrangian particle dispersion models, the 
McNider particle dispersion model and the NPS particle 
dispersion model now under development. The fractal dimension 
and Shannon entropy uncover weaknesses in both models which 
are not obvious with standard geophysical measures. They also 
reveal similarities and differences between the atmospheric 
models and simple chaos systems. Combined, these chaos 
measures may lend detailed insight into the behavior of 


Lagrangian Monte Carlo dispersion models in general. 
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THESIS DISCLAIMER 


The reader is cautioned that computer programs developed 
in this research may not have been exercised for all cases of 
interest. While every effort has been made, within the time 
available, to ensure that the programs are free of 
computational and logical errors, they cannot be considered 
validated. Any application of these programs without 


a@aaicional verification is at the risk of the user. 
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I. INTRODUCTION 


A. OVERVIEW 

Unlike the relatively simple linear relationships studied 
in basic physics, natural phenomena often exhibit complex 
nonlinear behavior. One familiar consequence is the inability 
of meteorologists to make accurate long-term weather forecasts 
(Devaney, 1990). Like turbulence in other systems, atmospheric 
weather systems seem to follow patterns which are chaotic and 
unpredictable in the long term. The behavior of a particle 
released in such a system is deterministic rather than 
stochastic. However, small differences in initial conditions 
result in particle paths which diverge exponentially with 
time. Thus, since the resolution and accuracy of measurements 
of the initial conditions are always limited, forecasts will 
diverge from reality exponentially with time such that long 
term predictability becomes impossible (Gleich, 1987). 

In fluids, predictability relates to particle diffusion. 
Although the relationship is not simple, less predictable 
systems tend to more diffusive. 

A good Lagrangian particle model which can accurately 
simulate diffusion over a wide range of atmospheric conditions 
is of value in predicting toxic dispersion from various 


sources, such as biochemical warfare agents, Navy LNG tanker 


leaks, nuclear reactor and weapons plumes, industrial chemical 
plants, large rocket and missile launch emissions, aborts, and 
static test firings, and a variety of planned and emergency 
releases from tanks and hoses involved in storage, transfer, 
and transport of dispersive volatile toxins. Presently, 
predictions using existing particle models may vary greatly 
from reality when operating in certain atmospheric conditions. 
(Venkatram and Wyngaard, 1988). 

Heretofore, atmospheric particle dispersion models have 
generally been compared against turbulence data in terms of 
spectra and standard geophysical turbulence measures, such as 
the turbulent Kinetic energy (tke), its component velocity 
variances (0,7, o,?, and o,”), and the Brunt-Vaisala frequency 
(BVF), an atmospheric stability measure. The aggregate 
diffusion of many Lagrangian particles has also been compared 
with the results of Gaussian plume statistical theories as 
well as measurements of dispersing clouds (Venkatram and 
Wyngaard, 1988). However, complex but wholly predictable 
periodic behavior such as ocean tides may appear quite similar |. 
to truly chaotic or random behavior when using- such 
techniques. Another potential set of tools in analyzing 
particle dispersion model performance are the recently 
developed chaos metrics such as fractal dimension, D,, Shannon 
entropy, S, and the Lyapunov exponent, dA (Baker and Gollub, 


1990). 


B. OBJECTIVE 

This study examines the self-affine fractal dimension, D,, 
and compares it to the Shannon entropy, S, and the Lyapunov 
exponent, A. These metrics are first applied to the 1-D 
chaotic system known as the logistics difference relation and 
the chaotic 2-D Henon system. These simple well characterized 
systems are studied to determine similarities and differences 
between the above three chaos metrics. Two 3-D Monte Carlo 
Lagrangian scattering routines designed to Simulate 
atmospheric dispersion are then studied as representative of 


more complex real world systems. 


C. WHY USE CHAOS METRICS 

There are at least two good reasons to try chaos metrics 
on dispersion. One, chaos metrics may provide a means to 
discriminate between periodic wave behavior, chaos, and 
differing degrees of turbulence. Second, chaos metrics may 
provide some additional insight into the behavior of particle 
dispersion in 3-D scattering routines. 

Kamada and DeCaria (1992) have shown that though nocturnal 
periodic gravity waves have quite different dispersion 
characteristics from turbulence, the two cannot be 
distinguished readily by using standard atmospheric turbulence 
measures such as the Brunt-Vaisala frequency (BVF), vertical 


velocity variance (o,’), turbulent kinetic energy (tke), 


buoyancy length scale (1,), or spectral analysis using FFTs. 
The self-affine fractal dimension, D,, was shown to be the 
only facile wave/turbulence discriminant tested. Since the 
Shannon entropy and Lyapunov exponent are two other standard 
chaos measures, they might also be tested as possibly useful 
dispersion measures. 

In the 3-D Monte Carlo scattering routines which model 
atmospheric dispersion, such as the McNider Particle 
Dispersion Model, chaos metrics might quickly determine 
certain model characteristics for a given range of parameters. 
For example, in examining an expanding cloud of particles, the 
Lyapunov exponent would measure the divergence rate of the 
particles. The Shannon entropy would measure the evenness of 
distribution of particles within the expanding cloud. The 
self-affine fractal dimension could measure the total apparent 
distance a particle travels ina time, T, as a function of the 
time resolution, ¢€, used within T. This might give an 


indication of the range of scales over which mixing occurs. 


II. THEORY 


A. OVERVIEW 

There are several methods of measuring chaos’ and 
turbulence. Of interest for the present purposes are the self- 
affine fractal dimension, D,, the Shannon entropy, S, and the 
Lyapunov exponent, A. Some standard geophysical measures will 
also be used for comparison purposes in the 3-D case. 
Descriptions of each are given below and the logistics 
difference, Henon, and atmospheric particle diffusion systems 


to which they are applied are then described. 


B. FRACTAL DIMENSION (D,) 

The self-similar fractal dimension, D,;, can be described 
readily with the Cantor set (Devaney, 1990). To form this set, 
a straight line is drawn with the middle third removed, as in 
the second line of Figure 1. The two resulting straight lines, 
which are one-third the original line length, are then 
Similarly subdivided. The four new lines, all 1/9th the 
Original line length, are further subdivided, ad infinitum. 
From top to bottom in Figure 1, the resulting sets of lines 
are self-similar; that is, the manner in which the geometry is 
altered is repeated at all successive levels of resolution. 

The Koch snowflake is another geometrically self-similar 


figure (Devaney, 1990). It is constructed by removing the 


Figure 1, Construction of the Cantor Set. 


middle third of each side of an equilateral triangle, and 
replacing it with two pieces of equal length, creating a six- 
pointed star. This star has 12 sides, all of length 1/3. The 
middle third of each of the 12 sides is again removed, and 
replaced with two lines of length 1/9. Again, the process is 
repeated ad infinitum. 

Like the Cantor set, the Koch snowflake is also self- 
Similar {Figure 2}. The jagged sides appear geometrically 
Similar at increasing levels of resolution. 

The self-similar fractal dimension is defined by (Devaney, 


1990) as 


ln (Total length of pieces) (II-1) 
in (resolution) 


D3; = 
In the Cantor set example, the total length of segments of 


unit length at level n is 2", and the resolution level is 3", 


so that 


Zl 
pp ee eer 309). (II-2) 
ln 32 ideal 2 





Similarly for the Koch snowflake, there are 4 pieces for 


each level of n with a magnification of 3, so that 


; : 
D, = 2O 42 262 . (II-3) 


Alb ages 
With geometric figures, D, is unitless; the measure 
involves a length divided by a length. However this definition 


of fractal dimension cannot apply directly to a time series 


fe 


Construction of the Koch Snowflake. From Devaney 


Figure 2. 
(1990). 


trace, since it would involve the square root of the amplitude 


squared plus the time squared, i.e., 


ff 


ln Sy VAA? + Ne (II-4) 
Dz = 2 
: Io 

This definition must be adjusted when applied to time series 
data to ensure that the end result is independent of the 
arbitrary unit scaling between amplitude and time. 

There are other ways to characterize the fractal dimension 
of a system. For a more suitable times series measure, McHardy 


and Czerny (1987) redefined the length metric as 


L(e) = 2 IF(t+e) - F(t)| dt. (II-5) 


where F is the amplitude of the time series at time, t, ¢€ is 
the time increment, and T is the time window over which L is 
defined. 

This definition effectively avoids the units scaling 
problem. Since the inverse time in 1/é€ cancels the time units 
in the integral, L(é€) is only in units of amplitude. The 


fractal dimension is then defined as 


Ga inet (e) 


ITI-6 
dine ( 


D,=- 


McHardy and Czerny applied these definitions to the time 
variance of X-ray luminosity data from the Seyfert galaxy 


NG5506. Collins and Kastogi (1989) later applied McHardy and 


Czerny's definitions in analyzing gravity wave spectra from 
the atmospheric mid-troposphere. Recently, Kamada and DeCaria 
(1992) applied D, as a tool for discriminating waves from 
turbulence in nocturnal atmospheric boundary layer data. 

To actually compute this function, the integral is 


approximated with a numerical summation, so that 
St TF 
= Dr (t) 


Since at each resolution ét=e, the leading term on the right- 


hand side is always unity, so that the length is now 
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L(e) = } |F(t-e) - F(t)| . (II-8) 
0 


Thus L(€) is the total amplitude change over a time series of 
length, T, for a given time resolution, €. In this form it is 
quite clear that time is removed from the length 
determination, so that L(é€) does not depend on some arbitrary 


scaling between amplitude and time. 


C. SHANNON ENTROPY (S) 
1. Definition 
When a particle is first released, its initial 
location is Known and completely specified, so the information 
entropy (defined later) for its location is zero. Later, 
according to given equations of motion, its position diverges 


from the initial point. If the range of its possible positions 
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is partitioned into equal segments, then for a given time 
interval, T, its motion can be recorded in terms of occupancy 
time for each segment. Then the particle dispersion rate might 
be measured by the seeming degree of randomness of occupancy 
time or evenness of the state probability distribution over 
time period, T. Shannon or information entropy is defined by 
the state probability distribution, so entropy may be regarded 
aS a measure of dispersion rate for a time series of particle 
states. This can apply to the actual particle position, its 
velocity, or its phase velocity. Shannon entropy is defined 


Simply as (Baker and Gollub, 1990) 


N 
Ss = > y joi WS ELSE, 
i=l 
where 
S = system entropy, 
= number of permitted states, and 
N 
Pp; = probability of state i, such that ) p, = 1.0 
i=1 


Then, for N permitted states (or position intervals), 
the maximum possible entropy corresponds to equal occupancy 


time in each of the N states, so p, is 1/N for all i. That is, 


(II-10a) 


Expanding the summation results in 


(Gal 


a 1 a 
= ee Jhgiai + ca LTI-10e8 
Smax N N N ( 


= 
N 


lh 


4 
N 


and collapsing the common multiples gives 


jl cual 
g “wv (5, Inj) 
ae N N (II-10c) 


il 
| 
= 
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S 


max , 


which corresponds to total randomness, or a completely even 


particle distribution across all allowable states. Also, 
oe = 0 (II-11) 


which corresponds to all particles being in one state. 

For a randomly moving particle in 2-D space, the 
computation is similar. The 2-D space may be divided into, 
aun a 100 X 100 grid. If the particle is moving completely 
randomly, at time, t, it may be in any one of the grid boxes 


with equal probability. The entropy for this system is then 


4 a) _i_|+4 1 ln 1 ares eer 
10102 1002 BL, (0 La} 8} 


= 1902 | —+— mn 1 ' 
100? neo 


Co =a 
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whieh is thewsanewas 


S = ln #00 Seee. (II-12b) 
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2. Shannon Entropy for an N-D System 
For a 1-D system, the domain is simply partitioned 
into n intervals, and the probability computed for each 
interval. For a simple system such as recursion of the 


logistics difference equation, 
Se iy ial WP, Sa le ae <9 am Lm.) 


the probability for the ith interval, p,, is simply the number 
of times the interval has been occupied after n number of 
recursion steps, divided by the total number of steps. 

Note that the number of partitions should be 
appropriate for the total number of steps. Having too few 
intervals is equivalent to too low a resolution of the entropy 
and may result ina misleadingly low value of S. For instance, 
with only one interval, p, = 1, and S = 0. Ideally, the number 
of intervals for maximum resolution of the entropy of the 


AN 


system is e™, where is the positive Lyapunov exponent for 


the system and N is the number of steps (Baker and Gollub, 


1990, pp.126-129). Since \ is typically of order unity, e% 


can 
easily be a computationally intractable number. 

Again, for 3-D systems, assuming N equal segments 
along each axis, the number of partitions will be N’, so 
computation quickly becomes unwieldy for large N. In practice, 
N of order 107 to 10° has been used by some authors (Baker and 


Gollub, 1990, pg. 88) aS a compromise between entropy 


resolution and computational efficiency. In analyzing real 
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data, the number of partitions should be related in some 
direct fashion to the maximum resolution, ¢€;, of the measuring 
devices, and the total measurement time, T. The total number 
of partitions should probably not exceed T/€é€,, or otherwise 
even a completely random distribution would still result in 
some intervals unoccupied. 

3. Ln vs Log, 

Wolf maintains that the Shannon entropy should be 
computed with log, rather than the natural logarithm, log, 
(Wolf, 1986, pg. 276). With log,, the Shannon entropy relates 
directly to information in the form of binary bits, since one 
bit of information can be specified as being in only one of 
two states, e.g., true or false, on or off, one or zero. That 
is, the log, basis sets the entropy equal to the minimal 
length of binary code required to describe the state of the 
system. If all particles occupy only one state, turning the 
bit for that state to the "on" position is sufficient to 
specify the state of the system. If the particles are evenly 
distributed among all states, the length of binary code 
required to specify the system is equal to the number of 
states, which means that the system is completely random 
(Tribus and McIrvine, 1971). However, for the purposes of this 
study, the Shannon entropy can also be written in the 


following form: 
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N 
ESOT oboe reee (Ande ta) 


i=1 


The only difference in measuring entropy by the 
natural logarithm versus base two is the value imposed on the 
constant K. Most computations assume that K=1. As long as the 
method of computing S is similar when making comparisons of 


computed entropy, there need be no confusion. 


D. LYAPUNOV EXPONENT (A) 

For an expanding turbulent cloud of particles, one measure 
to describe the system is the particle divergence rate. Due to 
turbulence, the trajectories of two particles arbitrarily 
close will diverge with time. The divergence rate can be 
characterized by Lyapunov exponents. The Lyapunov exponent, Xi, 
is also a measure of the system's sensitivity to initial 
conditions. In an N-dimensional system, there will be N 
Lyapunov exponents; however, these do not ' necessarily 
correspond to coordinate axes. So in a Cartesian coordinate 
system, \,, A,, and A; are locally defined Lyapunov exponents 


which generally cannot be described as i,, A and j.,. 


yl 
Direction is adjusted for each point along the trajectory. In 
one-dimension, the Lyapunov exponent is defined by (Wolf, 


eo, Pp. 275) as 


ie 


N-1 
L = 1a le el (Il=1g 
n-o N 5% 


which can be described over a map domain as an integral, 


N 
fppin|e(a|ai. (I-16) 


0 


Xd 


Overall, the Lyapunov exponent can be viewed as a measure of 
the average local stretching rate of the particle trajectory, 
as indicated by the log of the length of the slope, |f’(i)|, 
weighted by ,the probability of encountering that slope. 
(Wolf, 1986) 

The Lyapunov exponent is related to the entropy, S, as 
well as the information loss rate. For instance, if a measured 
position is presently known to a precision of 16 bits, and 
AX = 2 bits per second, then the particle*s future trajectory 
cannot be predicted to any degree of precision beyond 8 
seconds (16 orbits). 

It should be noted that the Lyapunov exponent defines the 
average rate of loss of predictive power, and may vary locally 
along the orbit. 

The Lyapunov exponent is readily computed for one 
dimension, but in higher dimensions calculation becomes more 
complex. In three dimensions, the directions of trajectory 
divergence and contraction must be defined in terms of local 
tangents, which vary from point to point, so the calculation 


of the exponents must constantly adjust for each change in 
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direction. Wolf has developed algorithms for computing higher 
dimensional Lyapunov exponents which involve reorienting the 
major axis of an ellipse for each point, then renormalizing 
the function after every few points so that the unit ellipsoid 
does not overlap an attractor (Wolf, 1986). 

In 3-D, three Lyapunov exponents are required to classify 
the system. A negative exponent indicates a dissipative 
dimension. If all three exponents are negative, the system is 
dissipative, e.g., a pendulum settling down to a fixed point. 
If an exponent is zero, the system orbits about a fixed point. 
If one of the three exponents is positive, the system is 


chaotic; the orbital trajectories are diverging. 


E. GEOPHYSICAL MEASURES OF TURBULENCE 
1. Overview 
Geophysicists use several standard methods to describe 
turbulence. Discussed here are those methods used to examine 
the McNider and NPS dispersion models. 
2. Turbulent Kinetic Energy (tke) 


The mean turbulent kinetic energy, or TKE, is defined 


as 
tke = = (0% aceic.) oe (gil 75) 
where 
o, = standard deviation, u component of particle velocity, 
Oo, = standard deviation, v component of particle velocity, 
Oo, = standard deviation, w component of particle velocity. 
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3. Gradient Richardson Number (Ri) 


The gradient Richardson number is defined by 
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Ri = a -= (II-18) 
Za Z 


where the overbars Signify time averages, and 





earth*‘s gravitational acceleration, 

vertical position above the surface, 

potential temperature, i.e. the temperature 
normalized for adiabatic expansion with pressure, so 


DN YQ 
Hl wi I 


Ra 
9 = al Cp ; (LE=tey 
Po 


where R, is the gas constant for dry air, and C, is the 
specific heat at constant pressure. 

Ri is used in place of the Reynolds number as a 
dynamic indicator of turbulence when the atmosphere displays 
a non-neutral density profile. The numerator is proportional 
to the buoyant production or destruction rate of tke, 
depending on negative or positive Sign, respectively. The 
denominator is proportional to the shear generation rate of 
tke, which is nearly always positive. Thus, a positive 
Richardson number indicates a stable atmosphere (increasing 
potential temperature with height) and suppression of 
turbulence. Commonly, when the buoyancy destruction rate of 
tke exceeds 1/4 the shear production rate, turbulence is 


suppressed, i.e., when Ri > 1/4 (Stull, 1988, pg. 176). 
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4. Brunt-Vaisala Frequency (BVF) 
The Brunt-Vaisala frequency, BVF, iS a measure of 
static stability. BVF is defined by 


v (II-20) 


BVF ’ 
Oz 


Alla 


and in principal gives the highest gravity wave frequency 
which a fluid can support. It is undefined for negative 
temperature gradients, FCs: an unstable atmosphere. 
merbjan, 1989, pg. 35). 
5. Buoyancy Length (1,) 
The buoyancy length, 1,, is the standard deviation of 
the vertical velocity divided by the Brunt-Vaisala frequency, 


0 
1, = = (II-21) 
The buoyancy length is meant to be a measure of the dominant 
eddy scale. (Stull, 1988, pg. 310). 

6. Monin-Obukhov Length (L) 


The Monin-Obukhov length is given by 


L2- uz 6 (Tee 
gk w'6) 


where 


1) 


surface temperature flux , 
von Karman constant = 0.4 , 


3 
ee) 


a 
ll 


Oo 


and u, = 2 


oe (II-23) 


Mm 
Zo 


u. is the friction velocity, a measure of surface drag 
due to turbulent friction. The Businger function, y,, is a 
stability correction which is approximated with a rational 


function by (Kamada, 1992b) 


252 > 0(stabiey 


tly 


fo 
L 


Vv, = § isoo1secme 5.4151107 = 
: < O(unstable). 
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Dee... ae ne 
1 - 3.59002232 2 - 0.799168457(=) 


(II-24) 


In the convective boundary layer, L is proportional to 
the height at which the buoyant generation rate of tke matches 
the shear generation rate. L is the primary stability neaeene 
for the atmospheric surface layer and is also used in 
combination with the inversion height, z,, to characterize 
boundary layer stability. L > 0 indicates a stable boundary 
layer where vertical turbulence is suppressed by positive 
density gradients with height, z. L < 0 indicates an unstable 
surface layer where upward vertical motion is encouraged by 


negative density gradients. A stable atmosphere implies that 


20 


a vertically displaced air parcel will tend to return to its 


original height. (Businger, 1973). 


F. 1-D AND 2-D CHAOS EQUATIONS. LOGISTICS & HENON FUNCTIONS 
1. Overview 

The initial focus of this study is to compare and test 
some simple potential dispersion metrics. The fractal 
dimension, D,, the Shannon entropy, S, and the Lyapunov 
exponent, <A, are simple expressions which can monitor 
transitions between periodic and chaotic behavior, akin to the 
transition between laminar and turbulent behavior in a real 
fluid. Two such expressions, well characterized in the 


literature, are the logistics difference equation: 


De = WoC Les en) os Gi 2 5) 


and the Henon system: 


1 - ax. + y, : 


Xn+1 


(II-26) 
Vioeg, 2 
feemad and Tobochnik, 1988, pp.152-178; Baker and Gollub, 
1990). 
2. The Logistics Difference Equation 
Though simple and one dimensional, recursive iteration 
of the logistics difference equation results in chaotic 


behavior for certain parameters. For 0 < X% < 1andu< 1, xX, 


Benverges to 0. For 1 < # < 3, and the same range for Xp, x, 
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converges to p/4. For 3 < p < 3.6, x, fluctuates between 2° 
discrete points; while beyond uw = 3.6, the solution is nearly 
random, 1.e€., CHNaotie. [Bhigure..) 
3. Graphical Analysis of the Logistics Difference 

Equation 

The logistics difference equation is parabolic, so 
that for an initial x, (given as 0.04 with w=2.9 in Figure 4), 
GQrawing a vertical line to the parabola corresponds to x,. 
Then a horizontal line drawn to the inscribed straight line of 
Slope unity corresponds to recursing x,,, to a new value for 
x,- Reiterating this procedure marches the solutions to a 
final value of 0.655. 

This final value can be determined exactly from the 
original function. At steady state, x,,, = x,. Then for Figure 


4, where uw = 2.9, 
Mpa e 250 ee li iar (LI=25) 
and solving for x;,: 


xX, =0, «x, Stoo. (IT s26y 


n 

Note since the equation is quadratic, that there are 
two solutions. The x,=0O solution is not stable (as defined 
below). The general rule is: if the magnitude of the slope of 
the function in the region of the solution is greater than 
45°, the fixed point is an unstable solution; if the slope is 


less than 45°, the fixed point is a stable solution. 
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Figure 3. Logistics Difference Equation: 1-D Chaotic Motion. 
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When the function is chaotic, there are still only two 
solutions to the quadratic equation. However, both points are 
unstable: they both repel rather than attract the value of 


X,41- The two points can still be computed, resulting in 
oa Ope xX. = 0.74... (II-27) 


Another way of looking at this is by taking the derivative 


of x,,, with respect to x: 





=p(1-2x) . (II-28) 


Checking the chaotic system of Figure 5, when x, = 0.74, the 
slope is 


AX, +1 
ax 


yw 


= 3.9(1 ~2(0.74)) = -1.87 , (II-29) 





which is less than -1, and hence unstable. 

This concept is vital when visualizing what the Lyapunov 
exponent is measuring. Any time both solutions are at points 
where the slope is greater than 45° or less than -45°, the 
Lyapunov exponent is greater than zero (ln of the slope, which 
is greater than unity), and the system is chaotic. 

4. Bifurcation Diagrams 

A map of the possible values of x, for various yw is 
Called a bifurcation diagram [Figure 6]. With the bifurcation 
diagram, the complete behavior of the function can be 


determined. For example, for 0 < uw < 3, x, tends to only one 
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Figure 6. Logistics Difference Equation Bifurcation Diagram. 
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value; the repetitive recursion of x,,, to x, points to one 
attracting fixed point, corresponding to a dissipative system. 
At uw > 3, the bifurcation diagram splits into two fixed 
points. The solution of x,,, bounces back and forth between 
these two points. This is period doubling to period 2. Note 
for still larger values of m that there is another split to 
period 4, then period 8, followed by a rapid series of 
bifurcations culminating in chaos, where the points appear to 
be distributed over a wide range of x,. 

Another curiosity appears in the bifurcation diagram. 
There are numerous "windows" within the chaotic region. The 
most obvious is around uw = 3.82, where the period 4 and 8 
behavior seems to appear again. Increased resolution would 
show many more such windows. Also note that period widths get 
progressively narrower, and that for higher periods it is 
difficult to distinguish say period 64 from chaos. 

5. The Henon Equations 

As with the logistics difference equation, the Henon 
set has stable solutions for a range of parameters a and b 
which correspond to one fixed point. Other values for a and b 
result in two fixed points, or period doubling [Figure 7], and 
period 4 and higher. For still higher values of a and b, the 
result is chaotic behavior: the possible positions appear to 


be spread throughout a discrete range of x and y [Figure 8]. 
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The Henon bifurcation diagram resembles the logistics 
difference equation diagram, but has some obvious differences 
[Figure 9]. One, it does not show symmetric splitting, but 
rather is skewed downward from the centerline. Two, at the 
value, a = 1.08 (with b = 0.3), a window appears in the 
bifurcation map with new values for x,; the periodicity in 
this region occurs at points outside the previous range of x,. 
These new points seem to appear out of nowhere, and are not 


linked to previous points by bifurcation. 


G. 3-D ATMOSPHERIC PARTICLE DISPERSION MODELS 
1. Overview 

The second focus of this study is to examine chaos 
metrics as performance measures for 3-D Monte Carlo scattering 
routines. Of immediate interest is particle dispersion within 
the atmosphere. The McNider Dispersion Model is commonly used 
to estimate the transport and diffusion of atmospheric 
pollutants. The model estimates diffusion based on atmospheric 
flow predictions, and simulates the behavior of pollutant 
clouds by representing them as the ensemble trajectories of 
numerous Lagrangian point particles (Pielke, 1984). For this 
study, the entire range of possible atmospheric stabilities 
was studied, from -0.2 < 1/L 0.1. Therefore, to render the 
computations tractable, the flow predictions were obtained 
from boundary similarity theory developed by Sorbjan (1990) 


and Kamada (1992 a,b) rather than from a mesoscale windflow 
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model, since steady state mean flows could be assumed. The 
difference is not large, since most of the similarity 
algorithms are also used in the turbulence closure scheme for 
the windflow model. 

Several weaknesses were perceived in the McNider 
formulation. Therefore, a similar model was developed at NPS, 
and was also tested. It features a double Gaussian vertical 
velocity skewness scheme and damping of the particle 
oscillations within the boundary layer (Kamada 1992b). 

2. The McNider Dispersion Model 

For this study, mean flow components are neglected in 

order to focus on the turbulent fluctuations. The McNider 


Dispersion Model utilizes the form 


MCE Fea cea x (eb) Fe oe) AE, 
Vie toN Ge = yi bet vit) AE, (II-30) 
UiteeteN C)e=e7 (ber wl ty) At, 


where x, y, and z define the particle position in a Cartesian 
coordinate system, and the u, v, and w terms refer to the 
fluctuating turbulent velocity components, respectively. The 


velocity components are computed by (Pielke, 1984) 


Wie) sue Art) Rk CAG) tu (ee AC) ; 
Mie) = vier Atc)R (AG) + v(t — AE), (gine 331) 


Wie w(toneAee Re (At) + w(t = AG), 
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where R refers to the Lagrangian autocorrelation factors for 
each velocity component and is a function of At. The second 
terms on the right refer to random fluctuating components of 
the velocity, and are related to the flow field's turbulent 
kinetic energy. This random fluctuation of the velocity is 


constrained by the fluid's turbulent kinetic energy with 


o, = 6,1 - R2(At), 
0 = o/1 - R2(At), (11-3 
o1, O,/1 - R{(At), 


where the o terms refer to the standard deviations of the 


velocity components. These parameters are not directly 
computed from atmospheric windflow in McNider's model. Instead 


an empirical formula is used to deduce these values: 
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K, is the local exchange coefficient, described later in the 


NPS model discussion. To account for advection across vertical 
gradients of g,, 


Oo,, 

ao, = 6¢ = w(t-At) +a, , (II-35) 

and the horizontal components are determined by 

2 2 
iio ae BREE TS 

6, = 90, = | 2 | L (E26) 

SD, nO 

L 


The maximum wavelength in the vertical velocity spectra is 
given by 











a nl eee 
(0.55 + 0.382) 
§ 
= ies eae Ny eg OR Ag 
Persea - exp{ $2) - 0. 0003 exe $2] ppl) «1 EAR rg KG rg 
wae Fas 
ein—aea)) 
merez/L < 0, and 
I a Fit 2 eee eats aCe = >o . (II-37b) 


In the above, z, denotes the top of the planetary boundary 
layer. 


The parameter, Ri, is the gradient Richardson number, 


defined previously by equation II-19. The critical Richardson 


Bis, 


number, Ri,, beyond which turbulence is suppressed, is given 
by 
RI =O cede eee eae, (Ii=3e) 


where Az 1S measured in centimeters. The Lagrangian 


autocorrelation terms are determined by 











Z -At 
R, (At) = exp/ ee | f 
R,(At) = oxp/ =| (II-39) 
i 
ke =A EG 
ree ax t) ex =| 


T, is the Lagrangian integral or e-folding time scale for 
velocity autocorrelations, and is determined for each 


component from the turbulence spectra with 
1 2 (ray We 
T, = 0.2f6,A,/V ; (II-40) 
T,, = 0.2B,A,/V . 


where A, is the dominant wavelength for each component. 


Furthermore, the average velocity is defined by 


V = fut + v2 + We. (Tt >aae 


6, the ratio of Lagrangian to Eulerian time scales, is given 


by 
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Vv 
0, 
ea) ios (II-42) 
0, 
Vv 
Oy 


where 6 is restricted to B < 10. 


The horizontal components for peak wavelength are 


given ky 


(II-43) 


bin AlN 
Vv 
© 


The constant, A, varies as a function of stability in the 


unstable boundary layer, and is computed from 


a 


Zine a Zz vA Zz 
Oot (13 =) (1 = 15 =) 4" (0-55 0,38 —)., i a ba 
( 7! ( ae ( = Fal 


el 
3 


0.05(1 - 34) 3 (1 - 154)%5, ‘ig ee La (alee aL 
1g re L ie 


(II-44) 


Values of u’/(t-At) and v’(t-At) are determined randomly from 
a normal probability distribution with zero mean. However, the 
turbulence distribution of vertical velocity is skewed [Figure 
10], and the normal distribution is modified by changing the 


method of computing w: 
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aot + 2 - 7 , OF, 
wit Sao). = oe L (II-45) 
a (Gyr Ther ra ; 
o ie 
where 
a = -0.028 + 0.6 |P] , r. 
Hy = 0. 540g, ee 
and 
00+ 7 Sa 
0.68(1 - 15-23)6°°2- ee 
j= ( Z) L 
LZ zZ 
0.1 =a0m aC) cee feaser 6) 
(>) iG 
(II-47) 


The variables w* and w are random values obtained from a 
normal probability distribution with zero mean and standard 
deviation, o,. As noted by Pielke (1984, p. 179), this method 
is not entirely satisfactory, since it depends on the 
particular random number generator routine. 
3. The NPS Dispersion Model 

As mentioned above, L, the Monin-Obukhov length was 
proscribed for this study, while V, the mean windspeed at a 
height of two meters was set to 3m/sec and the surface 
vegetative canopy roughness length was assumed to be 0.05m. 
Other flow parameters required by the McNider and NPS particle 
models were computed from recently developed boundary layer 


parameterizations, given here and by Panofsky and Dutton 
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Figure 10. Probability density function of vertical velocity 
in the convective boundary layer. From Weil, Dispersion in 
the Convective Boundary Layer, Lectures on Air Pollution 
Modelling, Venkatram and Wyngaard, 1988. 


Sie, 


(1984), Sorbjan (1990) and Kamada (1992a). More detailed 
descriptions of the underlying theory will be published ina 
forthcoming article. Many of the algorithms listed here are 
actually a part of the mesoscale windflow simulation which is 
also used to drive the McNider model. From this windflow 
Simulation the McNider model obtains the following: the 
turbulent diffusivity, K,, gradient Richardson number, Ri, 


vertical windshear, the windspeed at height z, v,, and surface 


zf 


layer friction velocity, u.. 
us. is computed from L, using the intermediate Businger 


function, y,, as determined by equation II-24, and the square 


m7? 


root of the surface drag coefficient, 


Ci/2 . k 
a | (Il—4e9 
in 2) = 
(z 
so that 
= Le II-49 
u. = MAX(C{.°,0.0a8) ( ) 


Given L and u., one can compute the surface vertical 


temperature flux, 


3 
ar. _ UD (II-50) 
Kok 


This allows an estimate of the free convective boundary layer 


turbulent velocity scale (for L < 0), 
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1/3 
ol 
a a, 7a (II-51) 


This can be modified to include shear induced surface layer 


turbulence (forced convection) according to 


OF 4 - 2/8 - in 





1 - 150-29] -14 
Lw3 (II-52) 


oe 0.4 
(Kamada, 1992a). 
For z/L < -0.5, 1i.e., above the surface layer, the 
following forms were used. The vertical velocity variance was 
given by 


2/3 2/3 2/3 1/3 
ee 1.12) 2 : =) F wn{ Z| 2 Bias 0| 
Za vine ae Zz; 


a EE wi 


f 





(II-53) 


where R= 0.2 and D = 0.1 are ratios of the inversion/surface 
temperature flux and entrainment zone/boundary layer depth 
(Sorbjan, 1990). 

In a standard atmospheric windflow model with second 
order turbulence closure scheme, the turbulence kinetic energy 
(tke) is computed prognostically. If o,? is supplied as above, 
the horizontal variances follow by subtraction. Without such 
a model, o,’? and o,’ are obtained diagnostically as follows. 

The Andre (1978) third order turbulence closure used 


the following prognostic form for vertical velocity variance, 
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Wo eee 7 Oo, 2] 2 (II-54) 
oe = d + 2D Ween Cc ne ee 5S 


Assuming steady state and neglecting diffusion, this can be 
truncated to estimate the anisotropy, o,’/@ where @ is the tke. 
Mason and Thompson's (1987) large eddy simulation of neutrally 
stable flow showed that o,’/6 = 1/3 near the surface and 
increased with height. This occurs because the surface 
restricts vertical motion. Also, most boundary layer 
turbulence results from "surface no-slip" induced shear which 
creates mostly horizontal tke. The increase with height 
occurs because the shear drops rapidly, so o,,/ée gradually 
approaches the isotropic 2/3 value through pressure re- 
distribution, consistent with Mason and Thompson's results. 
Like the tke on which it depends, ¢€, the molecular 
dissipation rate of tke also decays gradually with height. For 


convective boundary layers, Sorbjan (1990) parameterizes € as, 


Sala 2) 
Se eae (II-55) 
4 
Ee = pew.2/Z, . (II-56) 
Putting the above together, the tke anisotropy ratio can be 


parameterized as, 
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on 2€) - € Bw'@!, 
SC g (I-57) 


e Se, E 





(Kamada, unpublished). The first term on the right hand side 
accounts for the height dependence under neutral stability, 
while buoyancy flux in the second term accounts for non- 
neutral stability. This form actually corresponds to the 
Lenschow et al. (1980) measurements better than the Andre 
model itself or derivatives thereof (Therry and Lacarrere, 


1980). 


So the horizontal turbulent velocity variances become 


(5/9)0,7( 2e/o,? - 1) ; and (II-58a) 


N 
l 


Osea (II-58b) 


This allows one to compute é, the turbulence kinetic energy 
(tke) as, 


6G=(0/+a07 + 4,7 )/2 : (II-59) 


The molecular tke dissipation length scale is parameterized as 


Pee= 2/(1/2 + 1.4€6'” ) (II-60) 


from which the turbulent mixing length is estimated as, 


ie=— MIN(Gl) o7/é, 2) , (II-61) 
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This finally yields the required turbulent diffusivity, 
K, = 0.441,6"" ; (II-62) 


(Kamada, 1992b). The buoyancy gradient for the unstable 


boundary layer is given by 


BOO/0Z = 0.60./2,((1-2/2,) 7 / (2/2) 7" - 
2R%(2/2,)°7/(1 - 2/2; + 4)™)  , 
(II-63) 
(Sorbjan, 1990). The above formulations were applied to the 
convective (unstable) boundary layer (z/L < -0.5) above the 
surface layer. 


For the unstable surface layer (L < 0, 
Diem, > =075), 


b. =4(1 297) (II-64) 
K, = kutz/o. (II-65) 
6 = (5.6K, /z.)’ : (II-66) 
Ri 2) (II-67) 
O. = Ol2 = 0eS52/ Gee ; 7 (II-68) 
2 = ’ (II-69) 
$,= (1-142/L )”" , (II-7op 
@. = (14.75% (-2/ ae (II-71) 
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and finally 
€,7 = ug,/Z, , where the subscript 0 refers (II-72) 
to the surface value. (Sorbjan, 1990). 
For the neutral to stable boundary layer and surface 


layer, where Le: 0, 


R 
le 


Zee O07 Prana (El —73 


wy) 
lt 


5) (6) (II-74) 


(Kamada, 1992b). So that, 


Camels -2(2 = 2/2,)~ ’ ite? 5) 
6, = ¢, , and (II-76) 
Orolo — 2/z,)~* ; (II-77) 


The local friction velocity and Monin-Obukhov lengths 


were characterized in Sorbjan (1990) as 


a = u.(l = 2/Z,)~ , and (II-78) 


ee (te /Z,) "7? (II-79) 


and the temperature flux at height, zZ, was given by 





weet wiO',(1 = 2/Z)5. . (II-80) 


The Brunt-Vaisala frequency at height, Zz, was given by 
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BVF = 4.3u.0.7(1 + 3.7Z/L) (1 = 2/2,)°""/2 


The tke dissipation rate at Z was parameterized as 
@. = 3.6(1 + 3.7Z/H) (lay 2 , and 
e= gue. 
The Richardson number was estimated as 
Zz 


Ri = ——————_ , 
5Z + ZL 


and the momentum diffusivity was given roughly by, 


KUsZ (iL Cana 


Lt aoe 


(II-81) 


(II-82) 


(II-83) 


(II-84) 


(II-85) 


(Sorbjan, 1990). Both the unstable surface layer and neutral 


to stable boundary layer cases used: 


O, = Ud, , 

i ge ’ 

oe, 

y, = 2Inf (1 + (1 -~i14zZ/)4) /2 7 9 ana 
6, = © + ©.( 1In(2/2,) eae , where 





0. = w'0',/u. 
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(II-86a) 


(II-86b) 


(II-86c) 


(II-87) 


(II-88) 


(II-89) 


The mean wind at height z and its shear were estimated to be 


Vets ( 11/2.) = Wo) / Kk , and (1I-90) 


ré 


eyoe = Uagleted.727/%)~*(1— 2z/z)~*/kz (II-91) 


(Panofsky and Dutton, 1984 and Sorbjan, 1989 ). 

The following details only the significant items 
distinguishing the NPS and McNider models. The NPS particle 
model utilized the above formulations for the velocity 
variances and tke rather than those from McNider. For the NPS 


particle model, the Lagrangian vertical timescale was also 


estimated differently according to 


To Navel le (L>0) , and (t= 92) 


w 


ie — 601 3 7/ wae (ual < ame) (II-93) 


Recognizing that the horizontal/vertical eddy aspect 
ratio decreases with increasing stability in the convective 


boundary layer, the horizontal Lagrangian integral timescales 


were estimated by 


i (2 - 40/L) T,  , and (II-94) 


u 


jb aaa a (II-95) 


The McNider algorithms for horizontal integral 


timescales were retained for stable cases. 
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Rather than use the McNider formulations, the vertical 
velocity skewness was simulated by converting every fifth 
updraft into a downdraft, while accelerating the updrafts and 
decelerating the downdrafts. So the updraft/downdraft 
probability ratio was set to 60%/40%, while the updrafts were 
50% faster than the downdrafts, in tune with atmospheric 
measurements. 

Non-stochastic buoyant forcing was also added to the 
NPS particle model for non-zero density gradients. For pure 


buoyancy driven motion, the force balance is just 


Ow 56 

— ———— ’ TI-96 

at 6 
where w is vertical velocity, and 66 is the change in 
potential temperature due to a displacement, 6zZ, away from the 
particle's neutral buoyancy height in a domain where 


d0/dz # 0. Here the neutral buoyancy height is taken to be 


the initial release height of the particle. Then, 


Ow=- 2 800t . (II-97) 


ae) 


If the particle's "memory" were perfect, then after a 


time, t, its buoyancy forced velocity would simply be 


| 


9 586 0t . (II-98) 


OS ct 
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However, to Simulate dilution by ambient fluid, the particle 
must gradually forget its neutral buoyancy level. Its 
mnemonic e-folding time will be T,, the integral eddy 
coherence time scale already used to compute the stochastic 
component. So at each time step the particle's memory dims by 


the LCactor, 





7 (II-99) 


In this case, after one time step, 


W, = -(9/8)66,6t_ , (II~100) 
as before. But for subsequent time steps, the solutions are 
Ww, = -(9/8) 6tC(Mé6, + 66,) , (ii~1019 
W; = -(9/8) St (M’Se, + M50, + 5@,) ; (II-102) 


or in general, 


m 
Ww, = -(g/9)6t © 60, Mt, (II-103) 
1 


This series grows quickly with m. However, also observe 


that 


Wi, = Mw, - (9/8) 68;,,6C . (117-104) 


t 


So only the last velocity need be retained. The buoyant 


forcing then added to the standard stochastic component 
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(adjusted for vertical velocity skewness) to obtain the total 
vertical velocity fluctuation. 
4. Random Number Generator Kernel 

The dispersion models utilize normally distributed 
random numbers to generate the fluctuating velocity. These 
fluctuations are scaled by the standard deviation of the 
velocity components to relate them to the turbulent kinetic 
energy of the flow field. The method used to generate random 
numbers in this study waS a variation of the congruence 
method, and is outlined in program McNider, subroutines NORNG 
and STRNUM (Appendix C). To check the distribution of the 
random number generation routine, the count distribution was 
plotted versus standard deviation [Figure 11]. It appears to 
have a generally Gaussian shape with perhaps slight skewness 
to the left of center for the 3600 iterations used in this 


test. 
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Figure 11. Random number distribution, Program McNider. 
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4.00 


III. ANALYSIS 


A. EQUIPMENT/SOFTWARE 

The logistics difference equation, Henon equations, and 
the McNider Dispersion model were all written in standard 
FORTRAN-77, and run on Intel 386-based personal computers 
using the PROFORT compiler. All graphs were produced with 


Golden Software's Grapher and Surfer programs. 


B. LOGISTICS DIFFERENCE AND HENON EQUATIONS 
1. Methodology in analyzing the Logistics Difference 
Equation 
Computer generated data from the logistics difference 
equation was produced using Program Feigenbm (Appendix 1). 
Since use of the logistics difference equation involves 
recursion, an iterative series rather than a time series is 


created, such that the fractal lengths were redefined using 


L(i) = 4 \F(i+e) -F(i)| di, (III-1) 


oy 


which is approximated numerically by 


N 
Li) = 9) eG Sera (III-2) 


i=k,k 


where 
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qy 
bs. 
MT 


value of the function at the ith iterate, 


iteration step, 


Mm. 
Ill 


and k = resolution in terms of number of iterations. 


Analogous to the Lagrangian particle models studied 
later, the resulting iterative series represents a single 
particle, shifting position, with each successive iteration 
corresponding to a single time step. To insure independence 
from initial values, the program was run for 3700 iterations 
and the first 100 points were discarded. N was set to 3600 
iterations for all plots, giving 3600 point data sets. 

Since recursion does not permit non-integer 
iterations, the values of k had to be integer divisors of 
3600. The 45 available values of k used are listed as an input 
data file in Program Feigenbm (Appendix 1). 

From a series of k values for L(k) for a given yu, a 
standard linear regression routine then determined D,. The 
Standard deviation of D,, was very small in chaotic 
Situations, and large for during period doubling. The reasons 
for this are discussed below. 

In computing D,, the last three values of L(k) and k 
were discarded. On most plots, the curve flattened for 
K > N/3. This was also noted by DeCaria (1992) in analyzing 
- time-series data. So in the parabolic logistics difference 
equation, one reason for discarding the highest k values is 


the rising number of fold-crossings with larger k. The 
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"pbarticle" trajectory is confined to the unit interval by the 
left-right symmetry or fold at mid-parabola. Therefore, the 
distance traversed during a fold crossing is not fully 
accounted as with a system without folds. The change in 
"apparent" traversed distance will be lessened as more fold 
crossings are involved. This results in decidedly non-linear 
slopes for log L(k) versus log k for k which are a sizeable 
fraction of the window. 
2. Analysis 

For values of wu corresponding to a fixed point, the 
plot of Log,L(k) vs Log,(k) is fairly flat, and the 
corresponding D, value is small [Figure 12]. For values of yu 
corresponding to period doubling (period 2), the plot shows 
large jumps from a baseline, corresponding to multiple 
harmonics of 2 [Figure 13]. The baseline corresponds to a zero 
length. For these plots, the zero was adjusted by adding one 
to avoid a baseline at - ~, a result of Log,,(0). This shift by 
one unit was important for later calculations; without it, D, 
> o for all periodic functions. 

A period 4 plot still shows regularity [Figure 14]. 
There appear to be three sets of overlapping Slope patterns: 
one, a baseline of slope zero; the other two with similar 
slopes but different amplitudes. Again, this can be ascribed 
to harmonics. The length is zero for every fourth data point, 


and so will remain zero for harmonics having k values of 4, 8, 
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Figure 13. Logistics Difference Equation. x,,, = 3.0x,(1-x,), 
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16, etc. The length is nonzero for every second data point, 
not including multiples of four. There is similarity in 
lengths for k values of 2, 6, 10, 18, etc. Likewise, there is 
Similarity in lengths for k values which are not multiples of 
feeee.9., 1, 3, 5, etc. 

At the onset of chaos (yu = 3.569946) as in Figure 15, 
all L(k) lift off the baseline; without periodicity, all 
lengths will remain non-zero. Thus, D, will increase while dp), 
decreases. Also noted is a sensitivity to initial conditions 
for identical values of mw, another characteristic of chaos. 
Figure 16, with x, = 0.04, is different from Figure 14 with x, 
= 0.1. However, Figure 15, with x, = 0.1, is virtually 
identical to Figure 17, with x, = 0.01. There may be some 
intrinsic pattern to this sensitivity to initial conditions. 
The fold crossing patterns might provide further insight. 
However, this question strays from the present focus. 

For full chaos, op, becomes relatively small, and the 
plot becomes almost linear, as shown in Figure 18. This plot 
highlights why the last three points were discarded in 
computing D,; the plot remains very linear sans large values 
of k. The main reason is the fold crossing patterns mentioned 
earlier. Additionally, not enough data points are retained in 
the length computation for large values of k to maintain 
Similarity. For k=1200 and N=3600, the computed length is the 


sum of only three distance measurements, which is not enough 
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to develop a good statistical mean length. 

A plot of D, versus pw shows many of the features noted 
in the previous bifurcation plot [Figure 19]. At uw * 3.0 D, 
jumps from zero (corresponding to a fixed point) to * 0.58, 
corresponding to period 2. At uw * 3.45 is seen a transition to 
period 4. Higher uw values display regions of periods 8 and 16. 
However, at this level of resolution period 32 or higher order 
bifurcations cannot be discerned. Within the chaos regime, 
about uw = 3.6, D, seems to oscillate, and shows several sharp 
peaks and dips. The apparent window at uw *~ 3.85 corresponds to 
a period 4 oscillation. 

A plot of entropy versus uw shows similar features 
[Figure 20]. Regions of period 2, 4, and 8 are readily 
apparent. As in the D, versus pw plot [Figure 19], differences 
between period 16 and above and chaos cannot be discerned at 
this level of resolution. The large window of low periodicity 
at uw * 3.85 is also seen, as are several other windows of 
periodicity at roughly pw = 3.63 and pw = 3.73. 

The Lyapunov exponent also shows the trends seen in 
the fractal dimension and entropy [Figure 21]. A dips sharply 
at uw * 3.25, 3.5, and 3.55, the centers for periods 2, 4, and 
8. Again, changes beyond period 16 are hard to resolve. In the 
period doubling region, A> = 0 indicates that the function is 
transitioning to the next bifurcation, a fixed point 


fissioning to two fixed points. For A} > 0, the function is 
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chaotic. There are also numerous "windows of periodicity" in 
the chaotic region which appear in the bifurcation, fractal 
dimension, and entropy plots. These windows are more readily 
identified by the Lyapunov exponent; if A drops below zero, 
the function is periodic. 

Figure 22 shows the strong correspondences between 
entropy, fractal dimension, and Lyapunov exponent. 
Periodicities and chaos are apparent with each metric, as are 
some differences. D, sometimes peaks where S and A) dip, as at 
wu * 3.725. The low values of XA and S indicate periodic 
regions, for which D, sometimes sharply increases. The 
bifurcation map [Figure 6] provides no clue as to the cause. 
However, if the fixed points are more widely separated than 
the average distance between successive points in the 
surrounding chaos regions, D, will become relatively large. 
This is less likely for real fluids because motion is then 
constrained by physical forces and conservation laws. However, 
it points to a possible schism between diffusion and 
dispersion rates even without mean flow. 

3. Henon Function Analysis 

Since the Henon equations are a 2-D extension of the 

logistics difference equation, Similar behavior is expected. 


With two dimensions, the length measurement is modified to 


67 


7 QO 
5.0 
apll 
4.0 
ore 
Ze 
ie 


0.0 


Chaos Metrics 


ee 


©) 


=ai0 


— 4.0 


one 


Figure 22 
Lyapunov 


NO 
CO 


ae oe. 5.4 5 oes 4.0 
578 


3 
: ee 
My + 3 
it S Vi 
i! i} 
t} d i 
vs 1 Noeeeemmnace es vt 
iy ; it 
i H ¥ 
if : 

fol i 

H 





2S S68 One 5 one bre 4.0 
ne 


. Fractal Dimension (D,), Shannon Entropy (S), and 
Exponent (A) for Logistics Difference Equation. 


Left axis \; right axis S and D,. 


68 


(ah SG Aen RE Te ee (III-3) 


In the period 2 case [Figure 7], a plot of Log,L(k) vs 
Log,,.k 1s very similar to the logistics difference equation for 
period 2 [Figure 23]. Again, for every k, evenly divisible by 
2, the length drops to zero. There are two sets of slope 
patterns: a baseline with slope zero, and an upper non-zero 
Slope (corresponding to non-even k-values). This upper slope 
appears to be quite straight, except for the largest k value. 

At the onset of chaos [Figure 24], again the baseline 
lifts, and the overall deviation in slope is’ smaller 
[Figure 25]. It is also noted that Log,L(k) only approaches 
zero when k is quite large. As noted before, for large k there 
are not enough distance measurements to adequately define a 
good statistical mean length; nor are there enough points to 
adequately define a periodicity. 

At full chaos [Figure 7], the plot of Log,L(k) vs 
Log,k is so linear that, in regions of full turbulence or 
chaos, D, can be defined with only two data points, i and o, 
corresponding to one small inner and one large outer scale of 


k. {Figure 26]. Then D, in this case is given by 
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This shows immediately that an increase in D, implies a 
relative increase in the apparent magnitude of the 
fluctuations discernible at higher resolutions, i.e., a shift 
in the amplitude spectrum toward smaller scales. This 
interpretation of D, will be of use in analyzing the 
atmospheric Lagrangian particle models. 

There appears to be a close correspondence between 
entropy S and D, for the Henon function [Figure 27]. Both show 
jumps corresponding to periods 2, 4, 8, and perhaps even for 
period 16. Both plots show corresponding changes when there is 
periodicity within the chaos region. 

A higher resolution look at the region, 
1.052 < a < 1.082, displays some interesting complexity 
[Figure 28]. Not only are there normal bifurcations, but at 
a * 1.062 a new branch appears with no obvious connection to 
previous points. This does not correspond to the fissioning 
process observed earlier, but rather to an entirely new set of 
solutions. At a * 1.080, points from this new branch break 


away and drift toward the original branch. 
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PMiciosrejLon nl Ootmec sa t-062, DD, Jumps in 
amplitude whenever the bifurcation map shows a window of 
periodicity [Figure 29]. D, based on x rather than the total 
distance r portrays much but not all the same information: 
jumps correspond to most of the same windows of periodicity 
[Figure 30]. The strong peaks in D, correspond to sharp dips 
in S; however, sharp dips in S do not necessarily correspond 
to sharp peaks in D, [Figure 31}. So a jump in average 
distance traversed between iterations always accompanies a 
drop in position randomness but not vice versa. This suggests 
that periodicity in the Henon system sometimes results from 
the sudden appearance of fixed points which are as closely 
spaced or more closely spaced than the mean distance between 
successive iterations in the surrounding chaos. 

Recognizing that there are two free parameters in the 
Henon equation (a and b), D, and S were mapped for values of 
O<as<si1.5 and0<b< 1.0 [Figures 32, 33]. The two 3-d maps 
are remarkably similar, and suggest that exp(D,) would be of 
the same order as S. The regimes of low periodicity are well 
defined. Where D, and S are zero, the Henon function has 
period 1, corresponding to one fixed attractor. The first 
jumps to periods 2 and 4 are well defined at this level of 
resolution. There are other steps of discernible width for 
high a and b values, on the left and back sides of the 


graphical "mountain". Another interesting feature is the 
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of Henon Function, with b = 0.3. 
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behavior for high values of b, where both D, and S show sudden 
increases in value. 

4. Conclusions drawn from Logistics Difference and Henon 

equation analysis 

The self-affine fractal dimension, D,, is a good 
discriminator of low frequencies in data, e.g., periods 2, 4, 
and 8 of a given length data set. However, at high 
frequencies, period 16 and above in the initial bifurcation, 
it is difficult to use D, to differentiate between periodicity 
and turbulence or chaos. The entropy, S, seems related to D,, 
but S measures the evenness of particle position distribution, 
while D, measures the changing apparent jaggedness of the 
function as the resolution varies. A sharp change in S always 
accompanies a sharp change in D,, but not always vice versa. 
Therefore, D, and S generally show the same trends within the 
chaos region, but not always. In such cases diffusion and 
dispersion paces are not equivalent. 

The Lyapunov exponent provides perhaps the most definitive 
information: for a given value of A, we know for certain 
whether the function is stable, periodic, or chaotic. 
Additionally, the Lyapunov exponent should be able to describe 
the degree of chaos: the greater the value of AX, the faster 
the orbital trajectory diverges, and therefore the more 
chaotic the function. Unfortunately, the Lyapunov exponent is 


not readily determined for multi-dimensional systems. 
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Analysis of these two simple systems, the logistics 
difference and Henon sets, shows that S, D,, and dA display 
correlated but not identical behavior, since they measure 
rather different things. These insights can now be applied in 


studying atmospheric Lagrangian particle models. 


C. MCNIDER PARTICLE DISPERSION MODEL ANALYSIS 
1. Methodology 

Four aspects of the McNider Particle Dispersion Model 
were studied: the divergences of vertical position, total 
position, velocity, and phase velocity. The position was 
measured in radial distance from the particle origin, which 
was arbitrarily set at x,y = 0, z = 100 meters. To simulate 
the real atmosphere, the boundary layer inversion height was 
varied linearly from 2,000 to 200 meters as 1/L, the inverse 
Monin-Obukhov length, was varied from -0.2 to 0.1. Mean 
velocities were set to zero and total particle velocity 
reflection was assumed at the ground surface and inversion. 

Hints as to model behavior are again provided by the 
bifurcation maps, which now depict the expansion and 
distribution of particle range with decreasing stability, as 
measured by 1/L. Model performance was checked against 
standard geophysical measures such as the Brunt Vaisala 
Frequency (BVF), turbulent kinetic energy (tke), and vertical 
velocity variance (o,”), as well as against the two readily 


calculated chaos measures, entropy, S, and the self-affine 
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fractal dimension, D,, for the range of 1/L values seen in the 
atmosphere. Using 100 partitions, the computed maximum 


possible entropy for all plots in the analysis was 


S, 


Danae Gee (III-5) 

For this study 1/L was set at a particular value, a 
Single particle was released and followed for 3,600 time 
increments of 1/6 seconds each (600 seconds total), then 1/L 
waS incremented and the analysis repeated, until the entire 
range was spanned using 1/L increments of 0.01. 

2. McNider No~Skew Routine 

For real unstable atmospheres, the vertical velocity 
distribution is negatively skewed [Figure 10]. Updrafts have 
higher velocities and thus occupy less volume than downdrafts. 
However, for the initial analysis of the McNider model, the 
vertical velocity turbulence distribution was not skewed. This 
was done to check the model without the ad hoc method McNider 
developed to introduce skewness, and which as outlined in 
Pielke (1984, pp. 178-179) ‘appears to be incorrect. 

An r-‘ versus 1/L plot of the particle position shows 
that on the negative (unstable) side, D, and S follow roughly 
the same trend up to 1/L = 0 [Figure 34]. From the analysis of 
D, in the Henon system, this suggests that the ratio of large 
to small-scale vertical distances between points varies little 
in the unstable regime without overlaying skewness. However, 


D, again jumps suddenly when 1/L exceeds zero. Since vertical 


S2 


fluctuations tend to decay in general with increasing 
stability, one expects that S would decay and D, would tend 
upward steadily from left to right for real turbulence in the 
atmosphere. Again, the discontinuities near 1/L = 0 are 
probably due to discontinuities across the neutral transition 
in the McNider algorithms. The small scale fluctuations in §S 
with 1/L are probably due to the random nature of the particle 
paths. 

On the positive (stable) side of 1/L, both D, and §S 
show sudden dramatic step increases, followed by D, tending 
upward while S tends downward. Again, this highlights the fact 
that S and D, are not measuring the same thing. The sudden 
jumps in S in the McNider model seem antithetical to the fact 
that negative buoyancy tends to suppress vertical motion and 
render it oscillatory in real fluids. For growing positive 
values of 1/L, the decrease in S shows that the distribution 
of points within the domain becomes less uniform, while the 
rise in D, suggests that the distance between successive 
points decreases more slowly on the small scale than for the 
larger, longer time scale fluctuations. This is qualitatively 
consistent with measurements under increasingly stable 
conditions in real fluids (Stull, 1988). 

The bifurcation diagram {Figure 35] also shows that 
for 1/L < 0 the total range of positions steadily decreases 


with increasing stability, while for 1/L > 0 it hardly varies. 
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This is because the total position change includes both 
vertical and horizontal fluctuations, and horizontal 
fluctuations are much less dependent on stability. 

These plots reveal other weaknesses in the McNider 
model. Unrealistically sharp fluctuations seem to occur near 
the neutral stability transition (1/L = 0), probably due to 
the dichotomous nature of the McNider algorithms, one set for 
stable conditions, another for unstable. From equation II-18, 
it is also clear for isotropic turbulence that oa,’ = 2/3 tke. 
Stability suppresses vertical turbulence, while instability 
enhances it. So for 1/L < 0, o,*? should exceed 2/3 tke. Figure 
34 shows increasingly qualitatively incorrect ratios with 
increasing instability. o,? also drops to near zero at neutral 
stability (1/L = 0) due to very low mean windshear values in 
the mesoscale flow simulation. This may not be as significant 
a problem in the prognostic windflow models for which the 
McNider model was originally designed as it is in the 
Similarity-based flow model simulator used in this study. 
However, an artificial "dip" would probably still appear. 

The bifurcation map of the total 3-D position seems to 
show a fairly constant range for 1/L > 0, a jump at zero, and 
a varying range for 1/L < 0 [Figure 35]. The plot density is 
too heavy to discern actual distribution patterns for given 
values of 1/L, although the computer screen displays a Moire- 


like pattern. What can be extracted from this plot is the 
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general density of points and the upper and lower boundaries 
(in magnitude) of position. 

In contrast to Figure 34, a plot of only the vertical 
position change (Az) reveals that the pattern for D, is quite 
Similar to that for S [Figure 36]. The bifurcation map shows 
an opposite jump in the vertical position as 1/L crosses the 
zero-point [Figure 37]. 

The velocity distribution for various values of 1/L 
seems fairly constant, so that the S curve has only a slight 
upward trend [Figure 38], except at values near 1/L = 0. 
However, D, tends downward strongly from left to right, jumps 
at zero, then again tends downward at around 1/L = 0.7, where 
it begins increasing in value. This upward trend in D, is 
probably caused by increasing negative buoyancy. The 
bifurcation map [Figure 39] shows that the velocity range 
tends generally downward, with a jump near 1/L = 0, and a 
change in trend at 1/L = 0.7. 

The phase velocity plots are almost identical to the 
total position plots. This is because in calculation phase 
velocity the velocities are small in magnitude compared to the 
position values. 

3. McNider Model, Skew On 
A glance at Figure 42 shows a big problem with the 


2 


vertical velocity skewness algorithm. The tke and o,° values 


are much too large. An expected maximum o,? will be ~ 10m’/s’, 
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Figure 36. McNider Particle Dispersion Model. Test of 
vertical position, Az, with no skew. 
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Figure 38. McNider Particle Dispersion Model, Test of total 


velocity, no skew. 
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Figure 40. McNider Particle Dispersion Model. Test of Phase 


Velocity, with no skew. 
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while the plot shows ~ 100 m’/s’. It seems that the algorithm 
is incorrect for very unstable values of 1/L. However, the 
stable side (positive values of 1/L) shows well behaved, 
reasonable values for tke and a,?. 

Examining only the stable side (positive 1/L) of 
Figures 42-49, D, and S are surprisingly flat for velocity, 
but not for total 3-D position, vertical position change, or 
phase velocity. Also, the bifurcation plot of velocity range 
becomes quite small for positive 1/L (Figure 47). Oddly 
enough, the small values begin not at 1/L > 0, but 1/L > 0.02. 
This corresponds to a sharp drop in entropy and fractal 
dimension near this value of 1/L, suggesting that the skewness 
algorithm is responsible. However, re-examining Figures 36 and 
38 with "no-skew", entropy and D, drop sharply at 1/L = 0.02 
as well as at 1/L = 0, suggesting that something other than 
the skewness algorithm is responsible, perhaps the inherent 
dichotomy in equations I1I-44 or II-37a. 

Again, this sharp transition near neutral values of 
1/L does not seem to accurately reflect nature, but rather 
appears to be a discontinuity in the algorithm solutions. 

Figure 47 also shows that the velocity range is much 
too high for the unstable case; an expected velocity value is 
on the order of 10 m/s or less. Oddly, there are several 
windows where the velocity range drops to low values, and are 


roughly what is expected, i.e., for 1/L * -0.475. However, at 
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Figure 42. McNider Particle Dispersion Model. Test of total 
position, Ar, with skewness. 
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Figure 43. McNider Particle Dispersion Model Bifurcation 
Map. Test of total position, Ar, with skewness. 
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Figure 44. McNider Particle Dispersion Model. Test of 
vertical position, Az, with skewness. 
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Figure 45. McNider Particle Dispersion Model Bifurcation 
Map. Test of vertical position, Az, with skewness. 
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Figure 46. McNider Particle Dispersion Model. Test of total 
velocity, with skewness. 
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Figure 48. McNider Particle Dispersion Model. Test of phase 
velocity, with skewness. 
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Figure 49. McNider Particle Dispersion Model Bifurcation 
Map. Test of phase velocity, with skewness. 


LOS 


more negative values of 1/L the o,?/tke ratio indicates that 
nearly all the turbulence is in the vertical component which 
is not seen in the atmosphere nor is physically likely (Stull, 
1988). 

4. NPS Model, Skew On 

The NPS model employed the double Gaussian skewness 
scheme as described in the theory section. It avoids the 
problems seen with the McNider model skewness scheme regarding 
the inappropriate high values for velocity, tke, and oa,” for 
unstable values of 1/L. Also, the tke and a,” better reflect 
the behavior of real fluids: oa,?/tke exceeds 2/3 for unstable 
values of 1/L, but not for stable values. The _ sharp 
discontinuity in transition from unstable to stable 1/L is 
reduced somewhat. 

Curiously, there is a distinct pattern to the strong 
fluctuations in D,, S, tke, and a,’ for unstable 1/L but not 
for stable 1/L values. This was also observed in the McNider 
skewness routine. The cause has not been determined; it could 
be due to limitations in the random number generator or 
perhaps an undiscovered program error. This is most obvious in 
the bifurcation map of vertical position range, Az [Figure 
53], where there are sharp jumps. 

The phase velocity plots are similar to the total 
position plots in all cases examined. The reason is readily 


determined: the phase velocity is dominated by the total 
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position, which varies from near zero to 1000 meters, while 
the velocity varies from near zero to roughly 3 meters/second. 

As in the McNider plots there is little indication 
that the entropy is affected by a rising BVF, while D, seems 
to rise as BVF rises, most notably in the vertical position 
plots. The BVFs indicate a long period compared with the 
1/6 second timesteps. The entropy analysis will not detect a 
particle riding on a low-frequency wave with long period. In 
such a case the particle distribution over time may appear to 
be roughly uniform between a minimum and maximum value. This 
constraint is similar to the resolution requirements for 
Beropy computation noted earlier. If the At in the model were 
increased from 1/6 second to say 10 seconds, while still 
retaining 3600 iterations, the entropy might also show a 


variance with the BVF. 
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Figure 50. NPS Particle Dispersion Model. Test of total 
position, Ar, with skewness. 


108 


On 


1000.0 1000.0 


SOONG 800.0 
OD 
= 
6 600.0 600.0 
ia, 
o 
O 
iS 
2 400.0 400.0 
ae 

ZO ORG 200.0 


Cag 





Figure 51. NPS Particle Dispersion Model Bifurcation Map. 
Test of total position, Ar, with skewness. 
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NPS Particle Dispersion Model. Test of vertical 


position, Az, with skewness. 
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Figure 53. NPS Particle Dispersion Model Dispersion Map. 
Test of vertical position, Az, with skewness. 
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Figure 54. NPS Particle Dispersion Model. Test of velocity, 
with skewness. 
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Figure 55. 
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velocity, with skewness. 


114 


Gene, 





1000.0 
» 800.0 200.0 
QD 
= 
© 
eZ 
a 600.0 500.0 
O 
& 
O 
= 
» 400.0 400.0 
Y | ; 
© 
is 
ae 
200.0 200.0 
0.0 0.0 
=o =O = O10 0.1 


Figure 57. NPS Particle Dispersion Model Bifurcation Map. 
Test of phase velocity, with skewness. 


> 


Iv. CONCLUSIONS 

Chaos metrics such as the self-affine fractal dimension, 
D,, entropy, S, and Lyapunov exponent, A, are useful tools in 
the analysis of system characteristics and behavior. D, 
measures the jaggedness of a trajectory, or relative magnitude 
of small versus large scale fluctuations. The entropy measures 
the randomness of state distribution. The Lyapunov exponent, 
difficult to implement in higher dimensions, measures the 
divergence of trajectories, and leads to a predictability time 
scale. As in both the chaos and atmospheric systems, D, and S 
often correspond closely, but not always. The Henon equations 
also demonstrated that exp(D,) is often of the same order as 
Se 

These chaos metrics when used in conjunction with standard 
geophysical measures such as turbulent kinetic energy, tke, 
vertical velocity variance, o,?, and the Brunt-Vaisala 
Frequency, BVF, can be applied usefully to atmospheric 3-D 
particle dispersion models. The McNider model as listed in 
Pielke (1984) has an inherent weakness in the skewness of 
vertical velocity variance. This weakness leads to obviously 
incorrect values of tke, a,’?, positions, and velocities for 
unstable values of the inverse Monin-Obukhov length, 1/L. The 
McNider algorithms also display an inherent discontinuity as 


the inverse Monin-Obukhov length crosses the zero point, 
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reflecting a transition’from unstable to stable buoyancy. This 
sharp discontinuity does not correspond to the real 
atmosphere, which suggests that model performance suffers at 
small values of 1/L. 

The NPS particle dispersion model's measured performance 
seems to better reflect reality. it does not share the McNider 
model's weakness in the skewness of vertical velocity 
variance, and displays more realistic ratios of oa,” to tke. The 
NPS Particle Dispersion Model also reduces the discontinuity 
in transition from negative to positive values of 1/L. 
However, the NPS model still has a problem with unstable 
inverse Monin-Obukhov lengths (1/L < 0); it displays a pattern 
to the vertical position distribution not reflective of real 
fluid behavior. Since the NPS model was developed in response 
to the current study, there may still be problems in coding or 
in the random number generator. The NPS model requires further 
refinement to model particle behavior in a fluid 
realistically. 

Previous performance metrics for atmospheric particle 
models were only designed to measure aggregate particle 
behavior. The above chaos metrics parameters also offer some 
insight into the model behavior of individual particles. 
Results indicate that atmospheric dispersion models require 
further development at a fundamental level in replicating 


turbulent diffusion. For example, large fractal dimensions in 


7 


atmospheric Lagrangian particle models seem indicative of 
strongly stable, rather than unstable, conditions, contrary to 
what seems likely for real turbulence. This points to the 
neglect in such models of motions at scales other than the 
dominant scale determined by o,,,.- Thus, atmospheric 
Lagrangian particle models may simulate only single scale or 
at most a limited range of large scale diffusion processes. 
Fluctuations more in accord with real velocity spectra might 
display more realistic entropy and fractal behavior. As is, §S 
and D, behave oppositely at times, which indicates that 
diffusion and dispersion rates do not have equivalent meaning, 
even in the absence of mean windflow. 

Simple periodic behavior preliminary to chaos in classical 
bifurcatory systems is also apparently not entirely equivalent 
to laminar wave behavior prior to the onset of turbulence. For 
example, if the time resolution of the analysis is much 
shorter than the period, the Shannon entropy for an 
atmospheric Lagrangian particle model may be quite large for 
wave motion, but small for periodic behavior. This is because 
the number of possible states in wave motion is not limited to 
the amplitude extrema as in the periodic motion. 

One weakness of this study of Lagrangian particle model 
performance is the unavailability of real data for comparison 
purposes. As stated earlier, atmospheric data is normally 


obtained in the Eulerian rather than Lagrangian frame. This 
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suggests a need for development of Lagrangian-based sensing in 
ie ba Re| turbulence experiments. Position and velocity 
measurements of particle markers in a wide range of 
atmospheric conditions, both stable and unstable, as a 
function of time, with rapid response remote sensors, would be 
useful in developing better models of particle dispersion and 
diffusion. As yet, since the gathering of data from real 
fluids in the Lagrangian frame is not feasible, these study 
results suggest that more extended application of chaos 
metrics to Eulerian measurements of turbulence in real fluids 
is warranted in order to gauge the performance of Eulerian 
turbulence models. Such applications need not be restricted to 
small-scale phenomena, since mesoscale Rayleigh-Benard 
convection also displays transitions from laminar cellular to 
chaotic behavior (Agee et al., 1973). Model improvements based 
on such tests may then be extended to the Lagrangian frame. 
With regard to the Eulerian frame, both fractal dimension 
and Shannon entropy are simple calculations which can be 
performed real-time in situ with current sampling devices, and 
should be simpler to implement than FFTs, since there are no 
transform matrices. One area to consider when conducting 
measurements of the self-affine fractal dimension is ene size 
of both e€, the time increment, and T, the time window width 
over which the length is being defined. € has a minimum size 
dictated by the response time of the sensors, while the size 


of T is critical when looking for intermittent or low 
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frequency events, e.g., gravity waves governed by the Brunt- 
Vaisala Frequency. If €,, is too small, then the time 
increment will fall in the range where the velocity 
autocorrelation decay is not exponential but constrained to be 
parabolic due to continuity and the viscosity of real fluids. 
At the same time €,,, must be at least ~ 3 orders of magnitude 
smaller than T to establish statistical validity. Establishing 
the appropriate €,,, is crucial to distinguish waves from 
turbulence, and probably applies equally for distinguishing 
intermittent and coherent structures in general. (Kamada and 
DeCaria, 1992) 

Another area for further study is using the Lyapunov 
exponent in analyzing 3-D turbulence. Although difficult to 
compute in 3-D, A provides a definite test of chaos in the 
particle diffusion rate. A study of all three chaos metrics 
seems necessary to clarify the relationship of diffusion rate 
to dispersion rate in these models. 

Standard geophysical turbulence measures such as tke and 
o,? might also be applied to classical chaotic systems such as 
the logistics difference and Henon equations. The definitions 
of timescale and averaging time must first be established, 
e.g., One iteration equals one time unit. The corresponding 
tke would be zero for a fixed point, and definite magnitude 
changes would occur as the function bifurcates to periods 2, 


4, 8, and to chaos. 
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The utility of chaos metrics in analyzing the 3-D 
Monte Carlo based particle dispersion models also suggests 
possible utility for other laminar/turbulent phenomena. These 
might include plasmas, acoustics, and laser cavities. Solid 
state free electron gas systems may also find these chaos 
metrics useful in describing phonon and electron transport. 
Phonon resonance in crystal lattices also reminds one that 
particle behavior in lattice gases and cellular automata may 
be studied with such chaos metrics. These measures might also 
be useful in modeling gas or liquid phase complex chemical 
kinetics or radiation, which have long been simulated by Monte 


Carlo schemes. 
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APPENDIX A. LOGISTICS DIFFERENCE EQUATION ANALYSIS PROGRAM 


CH KH KKH KH KKH KKK KKK IKK KKK KKK KKK KK KKK KKK KKK KKK KKK KKKKKKEKKEKKKEKKKKKKEKEKK KK 


PROGRAM Feigenbm 


This program computes self-affine fractal dimension, Da, for the 
Feigenbaum recursion relation, 


x(i+t1l) = 4 lambda x(i)( 1 - x(i), 


as a function of lambda. Da is defined for the iteration series as: 


d Ln[{ L(i) ] T 
Da = - —————_ 7 where L(i) = (1/i)|  |F(itl) - F(i)|ai 
d Ln[{ i ] o 
N 
= SUM |F(i) - F(i-1)| 
i=k,k 


and k is an exact integer divisor of N. Here we choose 
N = 3,600 to give us at least three decades to use in computing Da. 
This also allows 48 exact integer divisors of N in the k array. So 


this gives us 48 points of i vs. L(i) in the linear regression for Da. 
HK KKK KKK KKK KK KKK KKK KK KK eK KK KKK KKK KKK KKK KKK KKK KKK KK KEKE KKK KK KK KKK KK KK 


qagagqaanqgqagagnqaganqananqaagaaganqgnqagngnaggnqaanaaa 


DOUBLE PRECISION lambda, lglngth(48, 100), lgk(48), y, sum, x 
DIMENSION y(0:3600), sum(48, 100), k(48), nyval(48), xy(2,48) 


& ans(16), Dalmbda(3,100) 
DATA k/l, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 6,7 18,920 ,) 24 
& 30, 36, 40, 45, 50, 60, 72, 75, 80, 90, 100, 103, 106, Gee 


& 116, 120, 144, 150, 180, 200, 225, 240, 300, 360, 4007345—47 
& 600, 720, 900, 1200, 5.18007 360cG, 

OPEN (1, FILE='feigenis.dat', RECL=255, STATUS='new', ERR=700) 
OPEN (2, FILE='feigenll.dat', RECL=255, STATUS='new', ERR=800) 


c 
y(0) = 0.0 
c Get Da for different values of lambda from 0 to l. 
m = 100 
DO 500 1= 1, m 
el = 1 
lambda = 0.01*el 
IF (lambda .gt. 0.9999) lambda = 0.9999 
x = 0.01 
c Create the iteration series for a given lambda 
N = 7200 
WRITE (1, *)' i va 
DO 200 i=1,N 
y(l) = 4.*lambda*x*(1l-x) 
x as ¥i( 1) 
WRITE ((17150) “ie 1) 
150 FORMAT (1x, 14, 3x, £8.5) 
200 CONTINUE 
c With iteration series, get lengths and self-affine fractal dimension 
WRITE (2, *)' Ln(e) in ere uc 
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c Get total length as sum of differences for different spacings of i 


da25 
300 


350 
400 
c 
ec 
Cc 


qgqaqgdaa 


500 


550 


700 


800 


c** 


qgqgqgqagaagqgagdgqagqgqggqogngoganganaanqaaan 


DO 400 3 = 1, 48 
sum(j,1) = O. 
DO 300 i = 3600 + k(j), 7200, k(j) 
sum(j,1) = sum(j,1) + abs( y(i) ~ y( 1i-k(j) ) ) 
IF (j} .-gt. 44) 


& WRITE(6,250) y(i),y(i-k(j)),sum(j,1),1,3,k(3) 
0 FORMAT ( 1x, 3£8.4,315) 
CONTINUE 
Tr(a «gt. 445 
& WRITE(6,*)'sum(',j,',1) = ',sum(j,1) 


lglngth(j,1) = alogl0O( sum(j,1l) ) 
lgk(j) = alogli0( float(k(j)) ) 
xy(1,3) = 1gk(j) 

xy (2,)) lglngth(j,1) 
WREDE@(2,350) Lgk(j)), lgingth(), 1) 
FORMAT (1x, £8.5, 3x, £8.5) 
CONTINUE 


Do linear regression to get Da from 1o0g10(k) vs. loglO [L(€)] values 
CALL STLNRG(n, nyval, xy, ans) 


-ans(2) contains Da 
ans(14) contains the standard deviation of Da 


Put lambda, Da, and oDa ina 3 by 1 array 


Dalmbda(1l,1) = -ans(2) 
Dalmbda(2,1) = -ans(14) 
Dalmbda(3,1) = lambda 
CONTINUE 


WRITE (3,550) ((Dalmbda(i,l), i = 1,3), 1 = 1,m) 
FORMAT (1x, 3(£9.5,2x) ) 

STOP 

WRITE (6,*)'error opening feigenis.dat' 

STOP 

WRITE (6,*)'error opening feigenll.dat' 

STOP : 

END 


kkk kkk KR KR RRR RRR RRR RRR RRR RRR RR KKK RK 
subroutine stlnrg (k, nyval, xy, ans) 


This subroutine computes the slope and other statistics for a 
linear regression with several y values for each x value or with 
one independent variable. 


argument use description 

k input Number of different x values 

nyval input Array of number of y values for each x 
value. Dimensioned k 

xy input array of x and y pairs (xl, yl, x2, y2, etc.) 
of dimension (2,k) 

ans output Array of results computed, dimensioned 16. 
ans 
1 Number 
*2 Slope 


3 Mean of y 
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4 Mean of x 

5 y-intercept 

6 Sum of y**2 

7 Sum of squares mean 

8 Sum of squares slope 

9 Residual 

10 standard deviation of x 

standard deviation of y 
12 standard deviation error 
13 standard deviation y-bar 

*14 standard deviation slope 
15 standard deviation y-intercept 
16 F-Ratio for slope 


qgqQagaqgagagqagadadnqnnggnadaa 
ra 
ay 
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DIMENSION nyval(1), xy(1), ans(16) 
sumx = 0.0 
sumy = 0 


-0 
sumx2 = 0.0 
ans(6) = 0.0 
sumxy = 0.0 
n= 0 
nx AL 


ny = 2 

DO 20 j = 1,k 

nyval(j) = 1 

m = nyval()j) 

n=n+mMm 

DO 10 i=i1, m 

sumx = sumx + xy(nx) 

sumx2 = sumx2 + xy(nx)*xy(nx) 
sumy = sumy + xy(ny) 

ans(6) = ans(6) + xy(ny)*xy(ny) 
sumxy = sumxy + xy(nx)*xy(ny) 


10 ny = ny + l 
nx = nx + nyval(j) + 1 
20 ny = ny + l 
en =n 
ans(l) = en 
Sl = ans(1)*sumx2 - sumx*sumx 
s2 = ans(1)*sumxy - sumx*sumy 


enl = en - 1.0 
en2 = enl - 1.0 


ans(2) = s2/sl 

ans(3) = sumy/en 

ans(4) = sumx/en 

ans(5) = ans(3) - ans(2)*ans(4) 
ans(7) = ans(3)*sumy 

ans(8) = ans(2)*s2/en. 

ans(9.) = ans(6) - ans(7) - ans(8) 


s4 = ans(9)/(en - 2.0) 

endsl = en/sl 

ans(10) = sqrt(1.0/(enl*ends1) ) 
ans(11) sqrt((ans(6) - ans(7))/enl) 
ans(12) sqrt(s4) 

al3 = s4/en 

ans(13) = sqrt(al3) 

al4 = s4*endsl 


ans(14) = sqrt(al4) 
ans(15) = sqrt(al3 + al4*ans(4)*ans(4) ) 
ans(16) = ans(8)/s4 
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RETURN 
END 
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APPENDIX B. HENON EQUATION ANALYSIS PROGRAM 


C KRRKKKKEKKEEKKEKKEEKEKKEEKEREKEEKKKEEKEKEEKKEEEKEEKEEEEKKEKEKREEKKREKKEEK 


qg00 


Qq00 


Q0Q0 


C3 
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PROGRAM CHAOS 


Written by Ray Kamada & Korey Jackson 


+ 


DOUBLE PRECISION xX(5000), ¥(5000), DA, A, B, 

DAERR, L(100), Bl, S 
INTEGER NMAX,1I1,JMAX,IFLAG, J1, NMAX1, COUNT, EPSIL(100) 
OPEN(2, FILE='LENGTH.DAT', STATUS='NEW', ERR=150) 


WRITE(2,*)' A B DA DAERR' 
OPEN(7, FILE='ENTROP.DAT', STATUS='NEW' ) 
WRITE(7,*)' A B Ss! 


Establish NMAX 


NMAX=3600+1 

NMAX1=NMAX+100 

JMAX=NMAX /2 

PRINT*, 'JMAX=' , JMAX 

IF (NMAX.GT.5000) GO TO 100 


Initialize all variables 


IFLAG=0 
CALL DIVISR(NMAX, JMAX, EPSIL, COUNT) 


Now set up a chaos generating routine... 


0 


A=1.0570D0 
B1=0.3D0 
I1=0 
IF(I1.GT.600)GOTO 55 
J1i=0 
A=A+0.0010D0 
B=Bl 
I1l=I1 + 1 
IF(J1.GT.30)GOTO 50 
J1i=J1+1 
B=B+0.05D0 
print*,'MAIN1 A, B,',A,B 
IFLAG=0 
CALL HENON (NMAX1,A,B,X,Y,IFLAG) 
print*,'henon complete',a,b,iflag 


Error trap for x, y out of range (diverging solutions) 


IF(IFLAG.EQ.1)GOTO 35 


Comment out above line & replace w/below line when B is fixed 


IF (IFLAG.EQ.1)GOTO 30 
if(iflag.eq.1.)goto 200 


Analyze the data 


CALL ANALYS (NMAX, JMAX,X,Y,A,B, DA, DAERR, COUNT, EPSIL) 
print*,'analysis complete', a, b, DA 

WRITE(2,45) A, B, DA, DAERR 

FORMAT (—F 7 .473x%,, FO.3,. 5%, FS. 57 ox eo. 5) 

print*, 'WRITE to file complete in analysis' 

CALL ENTROP(NMAX,X,Y,S) 
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ay, 
c 49 
50 
ce 54 
Cc 
55 
C 
C 


ey 
100 


146 
148 


150 
200 


WRITE(7,47) A, B, S 

FORMAT (F7.4,3X, F6.3,5X,F8.5) 
GOTO 35 

CONTINUE 
GOTO 30 


CONTINUE 
CLOSE (2) 
CLOSE (7) 
GOTO 200 


PRINT*, 'ERROR....NMAX Greater than 10000....Array size too small' 
GOTO 200 

PRINT*, 'ERROR...Can not open henon.dat file' 

GOTO 200 

PRINT*, 'ERROR...Can not open analys.dat file' 

GOTO 200 

PRINT*, 'ERROR...Can not open length.dat file' 

PRINT*, 'COMPLETED MAIN2' 

END 


C KKK KKKKKKKK KKK KKK KHER KEKE KKKE KEKE KE EKKKEKEKREKKKEKEKEEKKKEKKEKREKKKKEKKRKEEK 
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SUBROUTINE DIVISR(NMAX,JMAX,EPSIL, COUNT) 
INTEGER CHECK, CHECK1, COUNT, NMAX, JMAX, EPSIL(100) 
CHECK=0 

CHECK1=0 

COUNT=0 

NMAX=NMAX~-1 

CHECK=CHECK+1 

IF (CHECK.GT.JMAX) GOTO 1050 
CHECK1=NMAX / CHECK 

CHECK1=CHECK1*CHECK 

IF (CHECK1.NE.NMAX) GOTO 1000 
COUNT=COUNT+1 

EPSIL (COUNT) =CHECK 

PRINT*, 'EPSIL=',EPSIL(COUNT), COUNT 

GOTO 1000 

CONTINUE 

NMAX=NMAX+1 

RETURN 

END 


C KEK KKHEKKKEKKKEEKKEKEKKKEKEKKKEKE KEKE KEK KKEKEKEEKKEKREKKKEKEKEKEEKKKKEKEKREKKKKKEK 
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SUBROUTINE HENON(NMAX1,A,B,X,Y,IFLAG) 

DOUBLE PRECISION X(NMAX1), Y(NMAX1), A, B, CHECK 
INTEGER NMAX1, N, IFLAG, NMAX 
OPEN(4,FILE='HENON.DAT', STATUS='NEW', ERR=700) 


Initialize variables 


PRINT*, 'PASSED TO HENON OK' 
DO 300 N=1,NMAX1 


X(N) =0.0D0 
Y(N)=0.0D0 
CONTINUE 


MAIN COMPUTATION ROUTINE 


First initiate X(1), if desired. 


X(1)=-0.65d0 
Y(1)= 0.38d0 


WRITE(4,*)' N X Y' 
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DO 400 N=1, NMAX1 
X(N+1)=1.-A*X(N) *X(N)+Y¥(N) 
CHECK=X (N+1) 
IF (CHECK.GE.100.OR.CHECK.LE.-100.) THEN 
IFLAG=1 
PRINT*, 'PROGRAM OUT OF RANGE IN X' 
ENDIF 
Y (N+1)=B*X(N) 
CHECK=Y (N+1) 
IF (CHECK.GE.1000.OR.CHECK.LE.-1000.0) THEN 
IFLAG=1 
PRINT*, ‘PROGRAM OUT OF RANGE IN Y' 
ENDIF 
IF(IFLAG.EQ.1)GOTO 760 
IF(N.LT.101)GOTO400 
e WRITE (4,350) N, X(N), Y(N) 
350 FORMAT(I10,5X,E15.5,5X,E15.5) 
400 CONTINUE 


C 
C NOW SHIFT THE VALUES DOWN IN X(Y), Y(T) 
c 
NMAX=NMAX1-100 
DO 600 N=1, NMAX 
X(N)=X(N+100) 
Y (N)=Y(N+100) 
600 CONTINUE 
o 
C Return to the main program 
Cc 
PRINT*, 'COMPLETED HENON OK FOR',A,B 
GO TO 770 
c 


700 PRINT*, ‘Error in opening HENON.DAT file' 
C 740 £CLOSE(2) 
C7750 CONTINUE 
GOTO 770 
760 PRINT*,'X value out of range...A,B TOSSED', A, B 
IFLAG=1 
140 CLOSE(1) 
RETURN 
END 
C KHKKKKKKRKKKKKKKKK KKK KKK KKK KKK KKK KKK KK KKK KR KK KKK KKK KKK KKK KK KKK KKK 
SUBROUTINE ANALYS(NMAX,JMAX, X,Y,A,B, DA, DAERR, COUNT, 
+ EPSIL) 
DOUBLE PRECISION SCRAP1, SCRAP2, L(100), X(NMAX), SUM(100), 
+ Y(NMAX), EPS(100), DA, A, B, DAERR, SCRAP3, XY(200), 
+ ANS(16) 
INTEGER NMAX, I, J, CHECK, CHECK1,JMAX,COUNT, EPSIL(100), 
+ COUNT2, nyval(100), COUNT3 


VARIABLE LIST: 

NMAX=MAXIMUM NUMBER OF DATA POINTS 

COUNT=NUMBER OF EVEN DATA POINTS THAT WILL 
DIVIDE EVENLY INTO NMAX 

L(J)=LENGTH FOR A GIVEN EPSILON 

EPSIL(J)=WIDTH, OR AN INTEGER DIVISOR OF NMAX 

EPS(J)=EPSIL(J) IN INTEGER FORM 

DAERR=STANDARD DEVIATION OF DA (LOG FORM) 


qaqaaqaaaaagaana 


OPEN (3, FILE='ANALYS.DAT', STATUS='NEW', 
: ERR=9010) 
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INITITALIZE THE VARIABLES 


Q0Q0Q 


DO 4000 g=1, 100 
L(J)=0.0D0 
SUM(J)=0.0D0 

4000 CONTINUE 


c Now compute L(epsil) 


DO 7500 J=1, COUNT 
CHECK=NMAX-EPSIL(J) 
CHECK1=EPSIL(J) 

DO 7000 I=CHECK1, CHECK, CHECK1 
SCRAP1=X (I+1) -X (I-CHECK1+1) 
SCRAP2=Y (I+1)-Y (I-CHECK1+1) 
SCRAP1=SCRAP1*SCRAP1 
SCRAP2=SCRAP2*SCRAP2 
SCRAP2=SCRAP2+SCRAP1 
SCRAP2=DSQRT ( SCRAP2 ) 
SUM(J)=SUM(J)+SCRAP2 

7000 CONTINUE 
L(J)=SUM(J) 
7500 CONTINUE 
Cc 
C QUIK DATA PRINT 
DO 7550 J=1, COUNT 
L(J)=L(J)+1.0D0 
SCRAP1=L(J) 
L(J)= DLOG10(SCRAP1) 
SCRAP3=DBLE(EPSIL(J) ) 
EPS (J)=DLOG10(SCRAP3) 
WRITE(3,7545) EPS(J), L(J) 
7545 FORMAT(5X,F8.5,5X,F8.5) 
7550 CONTINUE 
Cc 
C Now do a linear regression to find DA 
c First put into a format to use the canned linear regression 
c subroutine. 
COUNT3=2*COUNT 
DO 7700 J=2, COUNT3, 2 


I=J/2. 

XY (J-1)=EPS(I) 

XY (J)=L(I) 
7700 CONTINUE 


Cc 
c COUNT2 IS THE COUNT WITH THE LAST 3 VALUES TOSSED FOR Da 
c computation 

COUNT2=COUNT=3 

CALL LINREG (COUNT2, nyval, XY, ANS) 

print*,'made it out of linreg' 

DA=ANS (2) 

DA=-DA 

DAERR=ANS (14) 

GOTO 9100 
9000 PRINT*, 'ERROR OPENING HENON.DAT FILE' 
9010 PRINT*,'ERROR OPENING ANALYS.DAT FILE' 
9020 PRINT*, 'ERROR OPENING LENGTH.DAT FILE' 
29100 print*,'made it to 9100° 

RETURN 

END 
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subroutine LINREG (k, nyval, xy, ans) 


This subroutine computes the slope and other statistics for a 
linear regression with several y values for each x value or with 
one independent variable. 


argument use description 


k input Number of different x values 


nyval input Array of number of y values for each x 
value. Dimensioned k 

xy input array of x and y pairs (xl, yl, x2, y2, etc.) 
of dimension (2,k) 


ans output Array of results computed, dimensioned 16. 


ans 
1 Number 
2 Slope 

3 Mean of y 

4 Mean of x 

5 y-intercept 

6 Sum of y**2 

7 Sum of squares mean 

8 Sum of squares slope 

9 Residual 

10 standard deviation of x 

11 standard deviation of y 

12 standard deviation error 

13 standard deviation y-bar 

*14 standard deviation slope 

15 standard deviation y-intercept 
16 F-Ratio for slope 


KHRKEKKHKEKKKKEK KEK KKK KKK KKK K KKK KKKKKKEKKKEKEKKKKKKRKKKEKKKRKRKKK 


DOUBLE PRECISION nyval(1), xy(1), ans(16), sl, s2, sumx, sumx2 


+ sumy, sumy2, sSumxy, enl, en2, s4, endsl, al13, al14 


sumx = 0.ODO 
sumy = 0.0DO 
sumx2 = 0.0DO 
ans(6) = 0.ODO 


sumxy = 0.0D0 


n= 0 

nx = 1 

ny = 

DO 20 5 = 1,k 

nyval(j) = 1 

m = nyval(j) 

n=n+tm 

DO 10 i=i1, m 

sumx = sumx + xy(nx) 

sumx2 = sumx2 + xy(nx) *xy(nx) 
sumy = sumy + xy(ny) 

ans(6) = ans(6) + xy(ny)*xy(ny) 


sumxy = sumxy + xy(nx)*xy(ny) 
ny =ny +1 

nx = nx + nyval(j) + 1 

ny =ny + 1 

en =n 

ans(l) = en 

sl = ans(1)*sumx2 - sumx*sumx 
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s2 = ans(l)*sumxy - sumx*sumy 


enl = en - 1.0 

en2 = enl - 1.0 

ans(2) = s2/sl 

ans(3) = sumy/en 

ans(4) = sumx/en 

ans(5) = ans(3) - ans(2)*ans(4) 
ans(7) = ans(3)*sumy 

ans(8) = ans(2)*s2/en 

ans(9) = ans(6) -—- ans(7) - ans(8) 


s4 = ans(9)/(en — 2.0) 

endsl = en/sl 

ans(10) = sqrt(1.0/(enl*ends1) ) 
ans(11) sqrt((ans(6) - ans(7))/enl) 
ans(12) sqrt(s4) 

al3 = s4/en 

ans(13) = sqrt(al3) 

al4 = s4*endsl 

ans(14) = sqrt(al4) 


ans(15) = sqrt(al3 + al4*ans(4)*ans(4) ) 
ans(16) = ans(8)/s4 

RETURN 

END 


C BERK KKEKREKKEKKEKKEKKEKKEKKEKKEKEKKEKKEKEKKEKKKEKKEKKEKKKKKEKKEKKEKKKEKKEKKKKKKEKEKEK 


C 23 Apr 92 reviSion (goes with O0<a<1.25 plot) 
Cc 
SUBROUTINE ENTROP(NMAX,X,Y,S) 
DOUBLE PRECISION X(NMAX), Y(NMAX), S, XMAX, XMIN, YMAX, YMIN, 
+ PROB(200,200), INVNMA, P, XRANGE, YRANGE 
INTEGER IX, IY, 1, J, NMAX 
C First find maximum and minimum values 
XMAX=0.0DO 
YMAX=0.0DO 
XMIN=0.0D0 
YMIN=0.0DO 
print*,'finding max and minimum values' 
DO 2000 I=1, NMAX 
IF(X(I).GT.XMAX) XMAX=X(I) 
IF(X(I).LT.XMIN) XMIN=X(TI) 
IF(Y(I).GT.YMAX) YMAX=Y(T) 
IF(Y(I).LT.YMIN) YMIN=Y(T) 
2000 CONTINUE 
C Now initialize the probability array 
print*, ‘initializing P array' 
DO 2105 I=1, 200 
DO 2100 J=1, 200 
PROB(I,J)=0.0DO 
2100 CONTINUE 
aes CONTINUE 
C Count the number of points in each unit area 
C First set up the scheme 
XRANGE=XMAX-XMIN 
XRANGE=2.00D2/XRANGE 
YRANGE=YMAX-YMIN 
YRANGE=2 .00D2/YRANGE 
DO 2110 I=1, NMAX 
IX=INT( (X(I)-XMIN) *XRANGE) +1 
IY=INT( (Y(I)-YMIN) *YRANGE)+1 
IF(IX.GT.200) IxX=200 
IF(IY.GT.200) IY=200 
PROB(IX, 1Y)=PROB(IX, I1Y)+1.0D0 
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2110 CONTINUE 
C Now normalize 
INVNMA=1.0D00/DBLE (NMAX) 
DO 2120 IxX=1, 200 
DO 221571Y=1, 6200 
PROB(IX, IY)=PROB(IX,1Y)*INVNMA 
ZS CONTINUE 
2120 CONTINUE 
C Finally, compute the entropy S 
S=0.0D0 
print*, ‘computing S' 
DO 2130 IX=1, 200 
DO 2125 IY=1, 200 
IF(P.ne.0d0)print*, 'P=',P 
P=PROB(IX,IY) 
IF(P.EQ.OD0O) GOTO 2125 
S=S-P*DLOG(P) 
2125 CONTINUE 
2130 CONTINUE 
C Now return to the main program to write S for this value of A 
RETURN 
END 
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APPENDIX C. PROGRAM MCNIDER 
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PROGRAM MCNIDE 


This program 
1) computes Monte Carlo particle velocity and position as a function 
of time step, t+*“*t, using the McNider algorithm, Kamada vertical 
velocity skewness variations, and Kamada mesoscale windflow and 
turbulence simulation model. 


2) computes the chaos metrics, Da (self-affine fractal dimension) 
for velocity, S (information entropy), and lambda (the Lyapunov 
exponent), as well as Ri, the atmospheric Richardson number, BVF 
(Brunt-Vaisala frequency) 1 (buoyancy length scale), ou,v,w, the 
velocity variances, TKE, the turbulence kinetic energy, Db 
(self-similar fractal dimension for position). 


RHEE KR KH KIKI KKK KKK KEKE KEKKKKKKKKRKKKRKRKKRKRKK 


Q #0gq gag aqgnNaAaanN0n0na a0 


DIMENSION nyval(45), k(45), ipts(100), bi(2,1800) 
DOUBLE PRECISION Linv, lglngth(45, 100), 1lgk(45), 
vert (0:3700), r(0:3700), r1(0:3700), r2(0:3700), 
el, p(1000), entropy(101), up, dn, lyap, pu, pv, pw, 
sum(45, 100), DavsL(8,100), pk, ans(16), xy(90) 
EQUIVALENCE (vert, r ) 
DATAsky I2,c, 4, 5S, 6, 8, 9, 10, 22, 25, 16, 18, 20, 24, 25, 
& 30, 36, 40, 45, 48, 50, 60, 72, 75, 80, 90, 100, 
& 120, 144, 150, 180, 200, 225, 240, 300, 360, 400, 
& 450, 600, 720, 900, 1200, 1800, 3600/ 
OPEN (1, FILE='McNideis.dat', RECL=255, STATUS='new', ERR=700) 
OPEN (2, FILE='McNidell.dat', RECL=255, STATUS='new', ERR=800) 
OPEN (3, FILE='DavsL.dat', RECL=255, STATUS='new', ERR=900) 
OPEN (4, FILE='SvsL.dat', RECL=255, STATUS='new', ERR=1000) 
OPEN (7, FILE='lyap.dat', RECL=255, STATUS='new', ERR=1100) 
OPEN (8, FILE='bifurct.dat', RECL=255, STATUS='new', ERR=1200) 
OPEN (9, FILE='initial.dat', RECL=255, STATUS='old', ERR=1300) 


QF Ly QM 


RKKKKKRKKKKK 


INITIALIZE 
REKKKKKRKRKK 
READ (9,30) 

30 FORMAT (///) 
READ (9,*) iis, ill, iDavsL, isvsL, ilyap, ibifurct, noskew, 
& iMCNid, iKamada, dt, ivert, iphase, iveloc, windspd2, z0O, 
keocabe, mine, iw, ix, Uptactor, adnfactor, dwkactor, zinitial, 
& icount, iposi 
WRITE(6,*)'iis, ill, iDavsL, isvsL, ilyap, ibifurct, noskew, ' 
WRITE(6,*)'iMcNid, ikKamada,dt,ivert,iphase,iveloc, windspd2,z0,' 
WRITE(6,*) ‘Start, m2, iw, ix,upfactor, dnfactor,dwfactor,zinitial' 
WRITE(6,*)'icount, iposi' 
WRITE(6,*) iis, ill, iDavsL, isvsL, ilyap, ibifurct, noskew, 
& iMcNid, ikKamada, dt, ivert, iphase, iveloc, windspd2, z0, 
S&S Start ,7m2, iw, ix, Upftactor, Adntactor, dwractor, zinitial, 
cVLCOUnG,, Lpos 1 

c Comments on input parameters: 

c lis ~ = 1 switch gives position and/or velocity file 


ao0a0 0 
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¢ tik = 1 switch gives Ln e€ vs Ln(L(e) file 

C iDavsL = 1 switch gives Da vs 1/L file 

c isvsL = 1 switch gives entropy, S, vs 1/L file 

c ilyap = 1 switch gives Lyapunov exponent vs 1/L file 

e ibifurct = 1 switch gives bifurcation disagram file 

c noskew = 1 switch turns off vertical velocity skewness 

c computation 

e iMcNid = 1 switch turns on McNider vertical velocity skewness 
c scheme 

c iKamada = 1 switch turns on Kamada vertical velocity skewness 
c scheme 

c dt = time step size in seconds 

c ivert = 1 switch gives vertical height file 

c iphase = 1 switch gives phase space vector file 

e iveloc = 1 switch gives velocity vector file 

(3) windspd2 = input proscribed wind speed at 2 meters 

a Zi = input proscribed fully turbulent boundary layer 

c height 

(e zO = input proscribed surface roughness height (~1/7) the 
c average height of surface roughness elements if 

C element 

@ density is between 10 and 50% (like most vegetation) 
c 

c Do some idiot proofing 

(4 


IF ( ivert .eq. 1) THEN 
0 


iphase = 
iveloc = 0 
iposi = 0 
ENDIF 


IF ( iphase .eq. 1) THEN 
ivert = 0 


iveloc = 0 
iposi = 0 
ENDIF 
IF ( iposi .eq. 1) THEN 
ivert = 0 
iphase = 0 
iveloc = 0 
ENDIF 
IF ( noskew .eq. 1) THEN 
iMcNid = 0 
iKamada= 0 
ENDIF 


IF ( iMcNid .eq. 1) THEN 
noskew = O 
iKamada= 0 
ENDIF 
IF ( iKamada.eq. 1) THEN 
noskew = 0 
iMcNid = 0 
ENDIF 
c Initialize variables 
ml = 1 
n = 3600 
en =n 
WRITE (4,*)' S liye: 
WRITE (7,*)' Lyap ‘7p 
WRITE (8,*) ' 1/L bifurcation pt' 
_ bo 500 1 = mil, m2 
c Let Linv, the inverse Obukhov length, 1/L = -kgw'0'/(O*ustar**3), be 
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Cc 
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g9qgaqa g 


varied in increments along the interval from -1 to l. 


em2 = m2 
Linv = start + 0.3d0*(1-1) /em2 
Zi = 2000 - 6000*(Linv + .2) 
r(0O) = 0. 

r1(0) = O. 


vert(0) = zinitial 


Z = zinitial 
= 0.001 

vdp = 0.002 
= 0.003 


kkk kk RRR 


GET iteration series for a given 1/L 
kkk kk kkk 


IF -{ 11s sne. 18"GOTO 150 
pu = 0.0 

pv = 0.001 

pw = -.01*alogl10(1) 
pw = -0.001 

time = 0.0 

rmax = 0.0 

r2max = 0.0 

rmin = 0.0 

r2min = 0.0 
sigU2sum = 0. 
SigV2sum = 0. 
sigW2sum = 0. 

zsum = QO. 

BVFsum = OQ. 

BLSsum = QO. 

sumdw = 0. 

tkebar = OQ. 

wbhuoy = OQ. 

sumw = OQ. 

Sumwbuoy = 0O. 
WRITE(6,*)'rmax2 =',rmax2 
rmax2 = Q. 

rmin2 = QO. 

iflag = 0 


DO 100 i= 1, n + 100 


CALL McNid (x, y, 2, udp, vdp, wdp, dt, 
iKamada, Linv, windspd2, deltaz, zi, 

sumwbuoy, BVF, BLS, sigmau2, sigmav2, 
pu, pv, pw, i, time, wk, ustar, iw, 
Sigmau, Sigmav, Sigmaw, Tlw, sumdw, 


DOR 


BVFsum = BVFsum + BVF 

BLSsum = BLSSum + BLS 

SsigU2sum = sigU2sum + sigmau**2 
SigV2sum = sigV2sum + sigmav**2 
sigW2sum = sigW2sum + sigmaw**2 
zsum = zsum + Zz 

bltime = 0.3*zi/(sigmaw*3700. ) 
eddytime = 0.01*Tlw 

IF (Linv .ge. 0.000) dt 
IF (Linv .ge. 0.000) dt 
hE etLiny .gt. 0.000) at 
IF (Linv .lt. -0.00) dt 
dt = 0.166666667 


eddytime 


fou uw tt 
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MIN (eddytime, 
MIN (dt, 0.25) 


noskew, iMcNid, 

zO, wbhuoy, sumw, Ri, 
Sigmaw2, tke, icount, 
i<,eupractor, iflag, 
dnfactor, dwfactor ) 


bltime) 


MIN (1.2*zi/(3700.0*wk), eddytime) 


time = time + dt 
cd WRITE (6,*) ‘time',time 
c get position vector, rl 
r1l(i) = sqrt( x*x + y*y + (z-zinitial)**2 ) 
c get velocity vector, r2 
r2(1) = sqrt( udp*dup + vdp*vdp + wdp*wdp ) 
c¢ get vertical position series 
IF ( ivert .eq. 1) r(i) = 2 - Zinitial 
c get position series 
IF ( iposi .eq. 1) r(i) = rl1(i) 
c get phase space position 
IF (iphase.eq. 1) r(i) = sqrt( rl(i)*rl(i) + r2(i)*r2(i) ) 
c get velocity series 
IF (iveloc .eq. 1) r(i) = r2(1) 
c get min/max range of the values 
IF ( i .gt. 100 .and. r(i) .gt. rmax ) rmax = r(i) 
IF ( i .gt. 100 .and. r(i) .1lt. rmin ) rmin = r(i) 
IF ( i .gt. 100 .and. r(i) .gt. rmax2) rmax2= r(i) 
IF i .gt. 100 .and. r(i) .1lt. rmin2) rmin2= r(i) 
e WRITE (1,50) 2,7. Vedi), 223) ee 
c50 FORMAT (1x, 14, 3(3x, £12.5) 
100 CONTINUE 
150 CONTINUE 
WRITE (6,*)'iteration series created for 1/L = ', Linv 
IF ( ibifurct .ne. 1 ) GOTO 240 
Cc KKK K KKK KKK K 


c GET BIFURCATION diagram, 
value 


c¢ put all the unique points in "bi" 
Ee KkeRKKRKRKRKKRKRKRKR KKK 
S 

m= 1 

bi(l,m) = Linv 

bi(2,m) = r(1999) 

Slop = 0.01*(rmax - rmin 
es IF (iveloc .eq. 1) slop 
e IF (ivert .eq. 1) slop 
c IF (iphase .eq. 1) slop 
(oy Check every r(i) 

DO 220 i = 2000, n+99 
c Check each r(i+l) against all 

DO 180 3 = 1, m 
up = bi(2,j) + slop 
dn = bi(2,j) - slop 
IF ( r(itl) .gt. dn 
S Since r(itl) lies beyond margi 
S unique point, so add it to the 
180 CONTINUE 

bi(2,m+1) = r(it+1) 

bi(l,m+l) = Linv 

m = m+l1 
220 CONTINUE 
oy IF ( iveloc .ne. 1 ) bi(2,m 

WRITE (8,230) (( bi(i,j), i 
230 FORMAT ( Ix, £9.5,05x, “£2.05 
240 CONTINUE 


WRITE (6,*)'bifurcation dia 
IF ( ilyap .ne. 1 ) GOTO 26 


C RRHKKKKKKKKKK 


c GET LYAPUNOV exponent for this ser 


QC KKKKKKKKKKKK 


using the last 1600 points. 


For each 1/L 


0.005*11. 
0.005*zi 
0.002*3000. 


iow tu 


unique points in bi 


-and. r(itl) .1t. up ) GOTO 220 
ns of all "bi", we assume it a new 
list of values in "bi". 


) = r(1999) 
= 1,2), j} = 1,m) 
) ; 


gram points found for 1/L 
0 


A] 
a 


Linv 


ies, 


E316 


qgqagqaqagaaaa 


n 
lyap = lim (1/n) sum ae i ex )} , where f£'(x ) su(l -2x ) 
n>0 i=1 i+l i+1 pio edl 


where we start from y(100) to avoid dependence on initial condition, 
x0. 


lyap = 0.0 
DO 250 i=zl,n 
e lyap = lyap + dlog( abs(mu*(1.0 - 2.*y(i+100) ) )) 
lyap = 0.0000000000 
250 CONTINUE 
e Normalize the lyapunov exponent by N - 100, the number of points 
S scanned. 
lyap = lyap/en 
WRITE (7,255) lyap, Linv 
255 FORMAT (1x, 2(£9.5,2x) ) 
260 CONTINUE 
IF ( isvsL .ne. 1 ) GOTO 290 
ey 
QC RRKKKKKKKKK N 
c GET ENTROPY value, S = - sum p log p , for this iteration series & 1/L 
QC KREKKKKKKKKK L=l i: e i 
c 
c Zero out the number of points in each bin 
DO 265 ii = 1, 100 
ipts(ii) = 0 
265 CONTINUE 
c Subdivide the interval into 100 sub-intervals, 
S 1i, and check all pts(i), from 101 - 3600. 
entropy(l) = 0.0 
DO 270 i = 101, n + 100 
ii = 99*( r(i) -— rmin2)/(rmax2 - rmin2) + 1.01 
a Make sure we don't create extra sub-intervals 
IF {ii .gt. 200) ii = 100 
cs Increment the ipts counter for sub-interval ii for every r(i) found 
c nT} Wana es 
ipts (11) 0 = "spes( ia) ter 
2/0 CONTINUE 
DO 285 ii = 1, 100 
e If no points are in sub-interval, ii, leave the entropy unchanged. 
IF ( ipts(ii) .eq. 0) THEN 
si = 0.0 
GOTO 280 
ENDIF 
e Otherwise compute probability of a point being in ii as the # of 
c points in ii divided by N, the total number of points scanned. 
pts = ipts(ii) 
prob = pts/en 
si = prob*alog (prob) 
e Since the probability is less than or equal to l, the log is 
S negative, so si is less than or equal to zero, so the entropy will 
c grow more positive, if the r(i) are scattered over more than one 
(2 sub-interval. 
280 entropy(l) = entropy(l) - si 
ZOD CONTINUE 
WRITE (4,255) entropy(l), Linv 
PeeG.ll. eq. 1) WRITE (2, *)* iggy Teje (DERN) 
WRITE (6,*) ‘entropy value found for 1/L = ', Linv 
290 CONTINUE 


IF ( iDavsL .ne. 1 ) GOTO 500 
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C KKKKKKK 


c GET Da. First get total length as sum of y(i) changes for varying 
Cc **x*kkk*k spacings of il. 
nn = 45 
DO 400 j = 1, nn 
sum(j,1l) = 1.0000 
DO 300 i = k(j), 3600, k(j) 
sum(j,l) = sum(j,1l) + abs( r(it100) - r(i-k(j)+100 )) 
300 CONTINUE 


c Get log of the total length 
lglngth(j,1) = dlog10( sum(j,l) ) 
C Get log of the spacing 
pk = k(j) 
lgk(j) = dlogl0( pk ) 
S Put these log values into an array and send to the linear regression 
Cc routine. 


xy(2*j-1) = lgk(j) 
xy(2*j) = lglngth(j,1) 
IF (ill .eq. 1) WRITE (2,350) lgk(j), lglngth(j,1) 


350 FORMAT (1x, £8.5, 3x, f£8.5) 
400 CONTINUE 
(a WRITE (6,*) ‘program passed statement 400' 
S 
e DO LINEAR REGRESSION to get Da from log10(k) vs. loglO [L(¢€)] values 
(o. 
S Use only spacings up to 1200 points for the linear regression. So 
Cc skip spacings 1800 and 3600. 
nn = nn - 2 
CALL STLNRG(nn, nyval, xy, ans) 
o 
c -ans(2) contains Da 
c ans(14) contains the standard deviation of Da 
c 
c Put 1/L, Da, oDa, BVF, BLS, ow?, tke in a7 by 1 array 
c 
ptsinv = 1./(n+100.) 
BVFbar = BVFsum*ptsinv 
BLSbar = BLSsum*ptsinv 
Ssigw2bar = sigw2sum*ptsinv 
tkebar = .5*(Sigu2sum + sigv2sum + sigw2sum) *ptsinv 
zbar = zsum*ptsinv 
DavsL(1,1l1) = -ans(2) 
DavsL(2,1) = ~-ans(14) 
DavsL(3,1) = Linv 
DavsL(4,1) = BVFbar 
DavsL(5,1) = BLSbar 
DavsL(6,1) = sigw2bar 
DavsL(7,1) = tkebar 
DavsL(8,1) = zbar 
WRITE (6,*)'Da value found for 1/L = ', Linv 
500 CONTINUE 
WRITE (3,*)' Da oDa Linv BVF BLS' 
& Pn ow2 tke zavg' 
WRITE (6,*)' Da oDa Linv BVF BLS' 
& ey ow2 tke zavg' 
WRITE (3,550) ((DavsL(i,l), 1 = 1,8), 1 = ml,m2) 
WRITE (6,550) ((DavsL(i,1), i = 1,8), 1 = mi,m2) 


550 FORMAT (1x, 5( £10.47, ix) 207-2 317 o) 
STOP 

700 WRITE (6,*)'error opening McNideis.dat' 
STOP 
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800 WRITE (6,*) ‘error opening McNidell.dat' 
STOP 

900 WRITE (6,*) ‘error opening DavsL.dat' 
STOP 

1000 WRITE (6,*) ‘error opening SvsL.dat' 
STOP 

1100 WRITE (6,*) ‘error opening Lyap.dat' 
STOP 

1200 WRITE (6,*)'error opening bifurct.dat' 
STOP 

1300 WRITE (6,*)'‘error opening initial.dat' 
STOP 
END 


RA KKKRKKKKRKKKKKKRKRKKKRKRKKKRKRKRKKRKRKRKRKRKRKRKRKKKKRKRKKKKRKRKKKRKRKRKKKRKKKKRKKKRKKKKRKKKRKRKKEK 


SUBROUTINE STLNRG (k, nyval, xy, ans) 


This subroutine computes the slope and other statistics for a 
linear regression with several y values for each x value or with 
one independent variable. a 


argument use description 

k input Number of different x values 

nyval input Array of number of y values for each x 
value. Dimensioned k 

XY input array of x and y pairs (xl, yl, x2, y2, etc.) 
of dimension (2,k) 

ans output Array of results computed, dimensioned 16. 
ans 
1 Number 
*2 Slope 
3 Mean of y 


4 Mean of x 

5 yrintercept 

6 Sum of y**2 

7 Sum of squares mean 
8 Sum of squares slope 
9 Residual 


*~*OONQQNQOQNONQaQNaNanNnaANNnNaNgaNnNNNnaNnNNnNNnNnNNgNnNnnNnNaNnNNANaANnA 


standard 
standard 
standard 
standard 
standard 
standard 
F-Ratio 


an 


deviation 
deviation 
deviation 
deviation 
deviation 
deviation 
for slope 


s(16) 


of x 

Sof y 

error 

y-bar 

slope 
y-intercept 


kak KKKKKKRKRKRKRKKKRKKRKKRKRKRKKRKRKRKRKRKRKKRKRKRKRKRKKRKRKRKKRKRKRKKKRKKRKKRKKRKRKRKRKRKRKRKRKKRKRKRKRKRKRKRKRKEK 


DIMENSION nyval(1) 
DOUBLE PRECISION xy(1), 


WRITE (6,*)' program entered linear regression routine' 


sumx = 0.0 
sumy = 0.0 
sumx2 = 0.0 
ans(6) = 0.0 
sumxy = 0.0 
n= 0 

nx = l 

ny = 2 

DO 20 3 = 1,k 


mo 


10 


20 


nyval(j) = 1 

m = nyval(j) 

n=n+m 

DO 10 i=i1, m 

sumx = Sumx + xy(nx) 

sumx2 = sumx2 + xy(nx) *xy(nx) 
sumy = sumy + xy(ny) 

ans(6) = ans(6) + xy(ny)*xy(ny) 
sumxy = sumxy + xy(nx)*xy(ny) 


ny =ny + 1 

nx = nx + nyval(j) + 1 

ny = ny + l 

en =n 

ans(1l) = en 

sl = ans(1)*sumx2 - sumx*sumx 
s2 = ans(1)*sumxy - sumx*sumy 


enl = en - 1.0 
en2 = enl - 1.0 


ans(2) = s2/sl 

ans(3) = sumy/en 

ans(4) = sumx/en 

ans(5) = ans(3) - ans(2)*ans(4) 
ans(7) = ans(3)*sumy 

ans(8) = ans(2)*s2/en 

ans(9) = ans(6) - ans(7) - ans(8) 


s4 = ans(9)/(en - 2.0) 

endsl = en/sl 

ans(10) = sqrt(1.0/(enl*ends1) ) 
ans(11) = sqrt((ans(6) - ans(7))/enl1) 
ans(12) = sqrt(s4) 

al3 = s4/en 

ans(13) = sqrt(al3) 

al4 = s4*endsl 


ans(14) = sqrt(al4) 

ans(15) = sqrt(al3 + al4*ans(4)*ans(4) ) 
ans(16) = ans(8)/s4 

RETURN 

END 


KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KEKE KKK KKK KKKEKKEKEKKEKKKKEKKEKKKEKKEK 


Cc 
Cc 
Cc 


SUBROUTINE McNid (x, y, Z, udp, vdp, wdp, dt, noskew, iMcNid, 


& IKamada, Linv, windspd2, deltaz, zi, z0O, wbhuoy, sumw, Ri, 

& sumwbuoy, BVF, BLS, sigmau2, sigmav2, sigmaw2, tke, icount, 

& pu, Pv, pw, i, time, wk, ustar, iw, ix, upfactor, iflag, 

& sigmaudp, sigmavdp, sigmawdp, Tlw, sumdw, dnfactor, dwfactor) 
this subroutine computes particle velocity and position for different 
values of Ri#, etc. 


KKK KKKKKKKKEKKEKKEHK EKER EKKKKKEKKEKKEKEKKKKEKEKEKKKKKKKKKKKKKKKKKKKKKKRKKKRKKKRKRKKEK 


REAL Km, L, lmdamu, lmdamv, limdamw 
DOUBLE PRECISION Linv, ruu, rvv, rww, pu, pv, pw 


c Get mesoscale and boundary layer flow parameters 
c Initialize variables 


RM — 


zi = 2000 — 6000*(Linv + .2) 


ZOVEZI. = 92/21 

zOvrL = z*Linv 
ziovrL = zi*Linv 
abszovrL = abs(zovrL) 


L = 1./Linv 

CALL MESOFLOW( z, zi, z0O, Linv, windspd2, ustar, Km, Ri, 
u, v, w, Theta2, dudz, dvdz, wk, D, R, wptpZ, BVF, iw, 
BLS, sigmaU2, sigmaV2, sigmaW2, tke, Wz, i, wk, ustar, 
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moO 


a0 


oO 


aan 


& 


bdthetdz,Thstar, Deltheta) 
P(e tteiw)WRELE (6, *) 'zOvrzi, ZOvrL, ZLOVEL ;ZO0VEZ1, ZOVrL, ZiovrL 


Get A, lambdas, ou,v,wdp, and buoyancy length scale, BLS 


& 


IF ( zovrL .lt. 0.0 ) THEN 
Sigmaudp = ustar*(12.0 + 0.5*zi/abs(L) ) **0.3333333 
IF(i.lt.iw)WRITE (6,*) ‘ou’ ‘=u*(12+.5zi/!L!)**1/3', 
Sigmaudp,ustar, zi,L 

IF ( abszovrL .le. l. 

A= 0.31*(1. — 3.*zovrL) **(-0.333333)*(1l. = 15.*zovrL) **0.25* 
(0.55 + 0.38*zovrL ) 

IF ( 0.1*abs(ziovrL) .gt. abszovrL .and. abszovrL .gt. l.) 
O-O05*(120—" 3. *20vrL)** (—07933333)-(lae— 15.*zovrL) **0.25 
MAX ( A, 0.07) 

MIN ( A, 0.18) 

lmdamu = 1.5*zi 

IF ( z .le. 0.1*zi .and. abs(zovrl) .gt. 1.) lmdamw = 5.9*z 

IF ( z .le. abs(L) ) Ilmdamw = z/( 0.55 + 0.38*zovrL ) 

eee (SO See le 2 ands oz 4 lt. zi<) 

lmdamw = 1.8*zi*( 1. - exp(-4.0*zovrzi) 
- 0.0003*exp( 8.O*zovrzi ) ) 

Sigmawdp = Km/ (A*lmdamw) 

IF(i.lt.iw) 

WRITE(6,*) ‘at z,ow*’ ‘=Km/(A*lmdamw) ',z,sigmawdp, Km,A, lmdamw 
ELSEIF ( zovrL .ge. 0.) THEN 

Sigmaudp = 2.3*ustar 

lmdamu = 0.7*sqrt(zovrzi)*zi 

es = 70.0 

IF (z.1t.205.0) es = 0.35*z 

lmdamw = MAX( z, 2.9*es) 

Ric = 0.25 

Rifactor = ((Ric - 1.25*Ri)/Ric)**0.58 

IF(i.lt.iw)WRITE(6,*) 'Rif=((Ric-Ri) /Ric)**.58',Rifactor,Ric,Ri 

shear = sqrt(dudz**2 + dvdz**2) 

Sigmawdp = 1.2*es*Rifactor*shear 

IF(i.1lt.iw)WRITE(6,*) ‘ow’ ‘=1.2esRi,shear', 

Sigmawdp, es, Rifactor, shear 

BLS = sigmawdp/BVF 

ENDIF 
Sigmavdp = sigmaudp 
IF (ikamada .eq. 1) THEN 


A 
A 
A 


Use values derived from Sorbjan and Kamada from mesoscale subroutine 


IF(i.1lt.iw)WRITE(6,*) 'ou?,0v?,ow?',sigmaU2,sigmaV2,sigmaW2 


Sigmaudp = sqrt(sigmaU2) 
Ssigmavdp = sqrt(sigmaV2) 
Sigmawdp = sqrt(sigmaW2) 


IF(i.1lt.iw)WRITE(6,*)'ikmda ou,v,w‘*',sigmaudp,sigmavdp,sigmawdp 


To get vertical integral length, see pps. 244-247 Lectures on Air 
Pollution 


Get 


Get 


aln = 0.41*z 

lmdamw = 0.27*sigmawdp/BVF 

lmdamw = 1./(1./aln + 1./1mdamw) 

IF(i.lt.iw)WRITE(6,*)‘'ikmda aln, lmdamw',aln, lmdamw 

ENDIF 
IF (noskew .eq. 1) tke = 0.5*(sigmaU2 + sigmaV2 + sigmaW2) 
lmdamv = lmdamu 
IF(i.1lt.iw)WRITE(6,*) '‘lmdamu, lmdamv, lmdamw', lmdamu, lmdamv, lmdamw 
mean velocity 
V = sqrt( u*u + v¥v + ww ) 
betas 
IF(i.lt.iw) WRITE (6,*)'u,v,w,Sigmaudp',u,v,w,Sigmaudp 
IF(i.lt.iw) WRITE (6,*) 'sigmavdp,sigmawdp',sigmavdp, sigmawdp 
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betau 
betav 
betaw 


0.6*V/sigmaudp 
0.6*V/sigmavdp 
0.6*V/sigmawdp 


d IF(i.lt.iw)WRITE(6,*)'at z,Bw=.6V/ow'‘',z,betaw,V,sigmawdp 


Get Lagrangian time scales 


an 


Tlu 
Til 
Tlw 


0.2*betau*lmdamu/V 
0.2*betav*lmdamv/V 
0.2*betaw*lmdamw/V 


IF(i.lt.iw)WRITE(6,*) 'betau, lmdamu,V',betau, 1lmdamu,V 


da IF(i.lt.iw)WRITE(6,*) 'Tlw=.2*8Bw*lmdamw/V',z,Tlw,betaw, lmdamw, V 


IF ( ikamada .eq. 1) THEN 

IF ( Linv .ge. 0. ) Tlw 

IF ( Linv .lt. O. ) THEN 
Tlw = 0.3*zi/wk 
AspectR = 
Tlu = 
Tlv = 
ENDIF 


AspectR*Tlw 
Tlu 


lmdamw/sigmawdp 


2.0 = 40.0*Linv 


d IF(i.1t.iw)WRITE(6,*) 'Tlw=lmdamw/ow**',Tlw,1lmdamw, sigmawdp 


ENDIF 

IF ( Tlu .lt. 0.02*dt ) Tlu 
IF ( Tlv .1t. 0.02*dt ) Tlv 
IF ( Tlw .lt. 0.02*dt ) Tlw 


Get Autocorrelation function 


an 


Ru 
Rv 
Rw 


exp( -dt/Tlu) 
exp( -dt/Tlv) 
exp( -dt/Tlw) 


aaa 


Get standard 
Sigmautp = 
Sigmavtp 
sigmawtp 
IF(i.1lt.iw) 


IF(i.1t.iw) 


qgaaa00 


ruu = 
rvv 
rww 


abs (pu/3.d0) 
abs(pv/3.d0) 
abs (pw/3.d0) 


IF(i.lt.iw)WRITE(6,*)'dt,TLu,TLv,TLw', 


0.02*dt 
0.02*dt 
0.02*dt 


dt, Tlu,Tlv,TLw 


IF(i.lt.iw)WRITE (6,*)'at z,Rw=exp(-dt/Tlw)',z,Rw,dt,Tlw 
IF(i.-lt.iw)WRITE (6,*) 'Ru,Rv,Rw, ikamada',Ru,Rv,Rw, ikamada 
deviations of velocity 

Sigmaudp*sqrt( (1. 
Sigmavdp*sqrt( (1. —- Rv**2) ) 
Sigmawdp*sqrt( (1. —- Rw**2) ) 


- Ru**2) ) 


& WRITE(6,*) ‘at z,aw’*‘=ow' ‘Vv (1-Rw2)',z,sSigmawtp, sigmawdp, Rw 


&SWRITE(6,*)'ou***,ov‘**,ow’**',sigmautp,sigmavtp, sigmawtp 
Get random normal fluctuation velocities, 


utp & vtp, at time, t-dt 


d IF(i-lt.iw) WRITE (6,*) 'ruu,rvv,rww',ruu,rvv,rww 


call NORNG( ruu, pu ) 
call NORNG( rvv, pv ) 
IF( pw .ge. 0.0 )wtp = 
utp = sigmautp*pu 
vtp = sigmavtp*pv 
IF ( noskew .ne. 1 
call NORNG( rww, 
wtp = 


pw ) 


d IF(i.lt.iw) WRITE 

d IF(i.lt.iw) WRITE 

da IF(i.lt.iw) WRITE (6,*)'w**’ 
GOTO 500 

250 CONTINUE 

da IF(i.lt.iw) WRITE 


d IF(i.lt.iw) WRITE 
if ( iMcNid .ne. 


Sigmawtp*pw + wbuoy 
(6,*)'pu,pv,pw',pu,pv, pw 
(6,*)‘utp,vtp,wtp',utp,vtp,wtp 


-and. Linv 


upfactor* (pw*(Sigmaw2 —- sigmawtp) + wbuoy) 


.lt. 0. ) GOTO 250 


ow’ ***pw',wtp, Sigmawtp, pw 


(6,*)'pu,pv,pw',pu, pv, pw 
(6,*) ‘utp,vtp,wtp.ucn, ven, veo 
1 ) GOTO 400 


c modify wtp for skewness according to McNider method 


iwp = 0 


142 


300 
c Get random normal fluctuation velocity, wtp, at time, t-dt 


Q a 


20 


d 


iwm = O 
CONTINUE 


call NORNG( rww, pw ) 

wtp = sigmawtp*pw 

IF ( wtp .ge. O. .and. iwp .eq. O ) wp = wtp 
IF ( wtp .ge. 0. ) iwp = 1 

IF ( wtp .le. O. .and. iwm .eq. O ) wm = wtp 
IF ( wtp .le. O. ) iwm = 1 

IF ( iwp .eq. O .or. iwm .eq. O ) GOTO 300 


Get S function 


& 


On = 0.2* zovrL**0.2 
O.1 + (0.6/( 0.68*(1.-15.*zovrL) **(-0.25) 


Bee ZOVEL «.<gt.. 0. )S 

ie zOvrL .le. 0. )$S 
ee on 2OVELG 1). ) 

IF(i.-lt.iw) WRITE (6,*)'iwp,iwm,S',iwp,iwm,S 


Get alpha and eta 


alpha = - 0.028 - 0.6*abs(S) 
eta = 0.54*abs(S) 


Get wtp(t-dt) !!! Pielke prints this on p. 178 as wdp(t-dt) CHECK THIS 
by looking at McNider paper!!!!! 


400 


qgqQgqqggaqagagagqaaa 


ga €00 94 


IF ( zovrL .le. O. ) wtp = alpha*wp + wm/alpha - eta 
IF ( zovrl .gt. O. ) wtp = wp/alpha + alpha*wm + eta 
IF(i.lt.iw) WRITE (6,*)'‘'alpha,eta,wtp',alpha,eta,wtp 
GOTO 500 

CONTINUE 

IF ( iKamada .ne. 1 .and. Linv .lt. 0. ) GOTO 500 


Try Kamada-Berkowicz double gaussian skewness approx instead of McNider. 
This approach assumes that mean absolute vertical fluctuation velocity 
is ~ 
w* to mixed forced/free convection. We let the up/downdraft volume 
ratio in the convective BL ~ 0.4/0.6, with mean updraft velocity 


0.55wk, vertically averaged over the BL depth, where wk generalizes 


1.080w and the mean downdraft velocity ~ 0.870w. We assume that the 


two gaussian velocity distributions are displaced from zero by this 
amount with 60% of the distribution in the downdraft volume and 40% 
in the updraft volume. We assume that departures from Gaussian are 
small for stable cases. (See pps. 199-201, Ch. 4 Lectures on Air 
Pollution Modeling) 


utp = pu*sigmaudp*sqrt(l - Ru**2) 
vtp = pv*sigmavdp*sqrt(l —- Rv**2) 
wbuoy = Rw*wbuoy - (9.801/Theta2)*Deltheta*dt 
sumwbuoy = Sumwbuoy + wbuoy 
call NORNG( rww, pw ) 
IF ( pw .ge. 0.0 ) ikount = ikount + 1 
PF (-1t.iwy WRITE (6,%*)° 1kount, icount’, 1.kount, icount 
IF ( ikount .eq. icount ) THEN 
pw = ~pw 
ikount = 0 
ENDIF 


The above sets an (icount + 1)/(2*icount) chance of updraft & an 


E. 
20 


& 


(icount - 1)/(2*icount) chance of a downdraft. 


g.-, if icount = 5, then updraft/downdraft chances are 60%/40%. 


CONTINUE 

IF( pw .ge. 0.0 )wtp 

IF( pw .le. 0.0 )wtp 

IF( pw .ge. 0.0 ) 
wtp = upfactor*pw*sigmawdp*sqrt(l - Rw**2) + wbuoy 

IF( pw .le. 0.0 ) 


upfactor*sigmawtp*pw + wbuoy 
dnfactor*sigmawtp*pw + whbuoy 


& wtp = dnfactor*pw*sigmawdp*sqrt(l - Rw**2) + whbuoy 


IF( pw .le. 0.0 ) wtp = upfactor*pw*sigmawdp + wbuoy 
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e IF( pw .le. 0.0 ) wtp = dnfactor*pw*sigmawdp + wbuoy 
sumw = sumw + wtp 
d IF(i.lt.iw)WRITE(6,*) 'wtp,pw',wtp, pw 
d IF(iflag .eq. iw) 
d & WRITE(6,*) ‘wbuoy,Rw, 60,sumw', wbuoy,Rw,Deltheta,sumw 
d IF(iflag .eq. ix) WRITE(6,*) ‘'sumwbuoy,sumw', sumwbuoy, sumw 
500 CONTINUE 
c 


c Get fluctuation velocities 


udp = udpold*Ru + utp 
vdp = vdpold*Rv + vtp 
wdp = wdpold*Rw + wtp 


IF (ikamada .eq. 1) THEN 


udp = udpold + utp - (1. - Ru)*udpold 
vdp = vdpold + vtp - (1. - Rv)*vdpold 
wdp = wdpold + wtp - (1. - Rw) *wdpold 
ENDIF 
d IF(i.lt.iw) 
d & WRITE (6,*)'at z,wdp=wdpold*Rwt+wtp',z,wdp,wdpold,Rw,wtp 


c Get position 


x = x + udp*dt 
y= y + vdp*dt 
= z+ wdp*dt 


Z 

IFo (Zz elt 10% 20 eee 
z= abs(z) + 10.*z0 
wdp = abs(wdp) 

ELSEIF ( z .gt. zi ) THEN 
Z= 2.0*zi - z 
z = MAX(z, 0.9*2zi) 
wdp = -abs(wdp) 


ENDIF 

udpold = udp 
vdpold = vdp 
wdpold = wdp 


iflag = iflag + l 
IF (iflag .gt. ix) THEN 


d WRITE (6,1110)time,Linv,Ri,udp,vdp,wdp,x,y,z 

dq WRITE (6,1120)Km, V, ustar, wpTpZ, sigmaW2, tke, Tlw, Rw, Wz 
iflag = 1 
ENDIF 


RREEKKEKKERKEEEKEKKEKKEKEKEKEKEKEKKEEKEKEKEKEKEEKEKEKEKREKKEKKEKEEKEEKEKEEKEKEKKKRKKRKKRKKRKKEKEEK 


c debug formats 

d1000 FORMAT (10x, 2g11.3) 

d1010 FORMAT (10x, 3f11.3) 

d1020 FORMAT (10x, 4f11.3) 

d1030 FORMAT (1x, i5 ) 

d1040 FORMAT (10x, 3g11.3, i5) 

a1050 FORMAT (10x, 4f11.3, i5) 

ad1055 FORMAT (10x, 5f11.3, iS) 

adi060 CONTINUE . 

G1110 FORMAT (1x,'’tm 1/D RISO ly wexeveze 7) fo. 2 toe oto. oe 

d1120 FORMAT (1x, 'Km V u* Q ow? e Tw Rw Wz ',£f6.1,5f6.2,£6.1,2f6.2) 
RETURN 
END 


RRKEKKKKEKKKEKEEKKKKEKRKEKKKRKEKKRKEKRKEKKKRKKEERKKEKKKKRRKRKKKRKRKRKRKRKRKRKKKRKRKRKKRKRKR KKK KE 


SUBROUTINE MESOFLOW ( z, zi, z0, Linv, windspd2, ustar, Km, Ri, 


& u, v, w, Theta2, dUtotdz, dvdz, wk, D, R, wpTpZ, BVF, iw, 
& BLS, sigmaU2, sigmaV2, sigmaW2, tke, Wz,i,wk, ustar, 
& bdthetdz, Thstar, Deltheta) 


c This subroutine computes the mesoscale flow parameters used as . 
c inputs to the standard mcnider particle calculations. It also computes 
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c some alternate values for the ou, v, w, and other quantities. 
REAL k, Ll, l1nzovzO, LovZi, Km, L, lepsilon, 1k 
DOUBLE PRECISION Linv 
itn, kpezae £ / 9.601, 0.4, 2.0, 0.00017 

a | IF(i.1lt.iw)WRITE(6,*)'in MESOFLOW i = ',i 

c Zero out some variables 


w 

< 

ny 
tou w 

O 


BLS 
c Initialize some values 

De='O.1 

R = 0.2 

zOinv = 1./2z0 

InzovzO = alog(z*zOinv) 

Ease ce ZOm 2 = Zoe ZO 

Zany = 1, /2 

Z2invy = 1. / 24 

L = 1./Linv 

GZio= 9.801421: 


zlovL = z1*Linv 
Theta2 = 290.0 
zOvL = z*Linv 
c Get psim & psih for surface layer. If stable, 
psim = -5.*zlovL 
psih = psim 
d IF(i.lt.iw)WRITE(6,*) 'psim=psih=-521/L',psim, psih, zlovL 
c If unstable, 
IF ( Linv .1t. 0.) THEN 
Phine—— (1. = 125. *zlovE)**(—-0. 25) 
d IF(i.1t.iw)WRITE(6,*) 'phim=(1-2821/L) **(-%)',phim, zlovL 
C Get psim by Kamada (unpub) 
Pre ZlOvien.Lt..0 sand. ZlovL .gt.--O-.O0l)) ZlovL =1=0.01 
d IF(i.lt.iw)WRITE(6,*)'just before psim = ' 
psim = (.0159651383 = 5.4151107*zlovn)7 
& (1. = 3.59002231*zlovL -0.799168457*zlovL**2) 
d IF(i.1lt.iw)WRITE(6,*) 'psim=(a-cz1/L)/(1-bz1/L-d(z1/L)?2',psim, zlovL 
c Also from Dyer and Bradley (1982) BLM 22, 3-19, we have for psi 
H 
d IF(i.1t.iw)WRITE(6,*)'just before psih = ' 
psih = 2.*alog( 0.5*(1. + rev, 1. - 14.*zlovL )) ) 
a IF(i.1t.iw)WRITE(6,*) 'psih=2Ln(45(1+V¥ (1-14z1/L)))',psih, zlovL 
ENDIF 
c Get windspeed drag coefficient and friction velocity 
d IF(i.1t.iw)WRITE(6,*)'just before sqrtCdn = ' 
sqrtCdn = 0.4/( alog(z1l*zOinv) - psim) 
d IF(i.1t.iw)WRITE(6,*) '¥ Cdn=.4/(Ln(z1/z0)-psim) ',sqrtCdn,z1,z0,psim 
ustar = MAX( sqrtCdn*windspd2, 0.01 
d IF(i.lt.iw)WRITE(6,*) 'u*=MAX(¥ Cdn*V2,0.01)',ustar,sqrtCdn,windspd2 
(os 
c For dispersion, get vertical profiles of relevant variables, using 
c similarity theory from (Sorbjan ,Structure of the Atmos BL, Ch. 4) and 
c (Panofsky and Dutton, Atmos Turbulence). First get oft used variables 
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zOvZi = z*ziinv 

zOovL = z20*Linv 

ZiovL = 21*Linv 

Lovzi = L*ziinv 

ustar2 = ustar*ustar 

ustar3 = ustar*ustar2 

wpTpO = -0.25508*Linv*ustar3*Theta2 

d IF(i.lt.iw)WRITE (6,%*) 'wO= =-.255u**30/L' ,wpTp0O,ustar,Theta2,L 


Tstar = -wpTp0O/ustar 
IF ( Linv 22... 0.) {2HEN 
Wstar = (gzi*wpTp0/Theta2) **.33333 


d IF(i.lt.iw)WRITE(6,*) 'w*=(gziw0/O) **1/3',wstar,gzi,wpTp0O,Theta2 
c Replace w* with wk for mixed forced/free convection (Kamada, 
e NPS-PH-92-007. Here, add only surface layer shear production. 
Ar = (1. - 150.0*Z200vL) **0.333333 - 1 
d IF(i.lt.iw)WRITE (6,*)'Ar= (1 - 150z0/L)**(1/3) - 1!' 
d IF(i-lt.iw)WRITE (6,*)Ar, zOovL 
d IF(i.lt.iw)WRITE(6,*)'just before Wk = ' 
Wk = wstar*(1. - ( 3.125 - 2.5*alog(Ar) )*Lovzi )**0.33333 
d IF(i.lt.iw)WRITE (6,*) 'Wk=w* (1-3.124-2.51n(Ar) )L/zi)**(1/3)' 
d IF(i.lt.iw)WRITE (6,*)wk,wstar,Ar,Lovzi 
Thstar = -wpTp0/wk 
c Decide what algos to use, based on 2/L. 


ENDIF 
IF ( zovL .1t. -0.5 ) icase = 1] 
IF ( zovL .1lt. 0.0 .and. zovL .ge. -0.5 ) icase = 2 
IF ( zovL .ge. 0.01 ) icase = 3 
GOTO ( 330, 400, 430 ) icase 
330 CONTINUE 
RREERKEKKKK 
c Get Km, & ou,v,w for forced/free convective BL. Ri and shear not 
needed. 
KKEKKRKEKKRK 
Cc Get potential temperature at height z above -L/2 
R = abs(R) 
ZONZi13 = ZOVZ1«*.35555 
ZOVZi23 = zovzil3*zovzil3 
tpl = 1.0 = zovzi 
tp4 = tpl + D 
tp413 = tp4**.33333 
tp423 = tp4**.66667 
Thstar = -wpTp0/wk 
tp23 = tpl**.666667 
R23 = R**.666667 
Oldtheta = Theta 
Theta = Theta2 + Thstar*( tp23/zovzil3 + R23*zovzi23/tp4l13 ) 
IF ( i.eq. 1) Oldtheta = Theta 
Deltheta = Theta - Oldtheta 
d IF(i.lt.iw)WRITE (6,*)'50,0,01d0', Deltheta, Theta, Oldtheta 
e temperature and buoyancy fluxes at height, z 
wpTpZ = wpTpO*(1. -—- 1.2*zovzi) 
B = g/Theta2 
BwpTpzZ = B*wpTpz 
Get vertical velocity variance at height z above -L/2 
SigmaW2 = wk**2*(1.1*zovzi23*tp23 + R23*zovzi23*tp423) 
IF(i.lt.iw)WRITE (6,%*) 'ow2=wk2(1.1(2/zi)**2/3(1-z/zi)**2/3 + ' 
IF(i.lt.iw)WRITE (6,*)'R**(2/3) (z/zi)**(2/3) (1-z/zi+D)**(1/3)' 
IF(i.lt.iw)WRITE (6,*)sigmaW2,wk,zovzi,R,D 
SigmaW = sqrt(sigmaW2) 
c Get TKE dissipation rate at height z above -L/2, using similarity 
tmp10 = Wk**3*ziinv 


Q 


ao 
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epsilonO= 0.75*tmp10 
phiEp = 0.75*tpl + R*zovzi 
epsilon = phiEp*tmp10 
d IF(i-lt.iw)WRITE (6,*)'e = @e*wk**3/zi',epsilon, phiEp,wk, zi 
(el Get ou2,ov?2 
WPSQOVE = 0.333*( (2.*epsilonO - epsilon)/epsilonO + 
BwpTpZ/epsilon ) 


d IF(i.lt.iw)WRITE(6,*) 'w‘2/e=(2€0-e)/3€0 +Bw'O'z/e' ,wpsqovE 
d IF(i.1lt.iw)WRITE(6,*)'e€0, Bw‘O‘z', epsilonO, BwpTpzZ 
SigmaU2 = 0.5555*sigmaW2*( 2./wpsqovE - l.) 
c SigmaU2 = MAX( sigmaU2, .33333 ) 
SigmaV2 = 0.8*sigmaU2 
d IF(i.lt.iw)WRITE (6,*)SigmaU2,wk, zovzi,sigmW2 
e The following tke parameterization is OK for low shear 
tke = 0.5*( sSigmaU2 + sigmaV2 + sigmaW2 ) 
d IF(i-lt.iw)WRITE (6,*)'e = ow2/2 + 30u2/2',tke,sigmaW2,sigmaU2 


tkel2inv = 1./sqrt(tke) 
tke32inv = tkel2inv**3 


c Get dissipation & mixing length scales & vertical diffusion 
c coefficient 
lepsilon = 1./( zinv + 1.4*epsilon*tke32inv ) 
d IF(i.lt.iw)WRITE(6,*)'le = 1/(1/z +1.4eve )',lepsilon,z,epsilon,tke 
lk = MIN ( z, 3.0*lepsilon*sigmaW2/tke ) 
d IF(i.lt.iw)WRITE (6,*)'1lk = MIN(3le*ow?2/tke, z)',1k,lepsilon,tke,z 
S The mixing length scale and diffusion coefficients are Kamada 
c variants 
Km = 0.44*1k/tkel2inv 
d IF(i.-lt.iw)WRITE (6,*)'Km = 0.44*lk*sqrt(tke)',Km,lk,tke 
c Get temperature and buoyancy gradients 
c dthetadz = .6*ThsStar*ziinv*(tp23/zovzi23**2 -2.*R23*zovzi23/tp423) 
dthetadz = .6*Thstar*ziinv* (tp23/zovzil3 =2.*R23*z20vZ123/tp423) 
Bdthetdz = B*dthetadz 
GOTO 460 


400 CONTINUE 
Kae KkK KKK KKK 
c get Km for unstable surface layer. Shear and Ri not needed for 
unstable. 
kkk KK KK KKK 
Km = k*ustar*z/phim 
tkeusl = (5.6*Km*zinv)**2 
Ri = zovL 
phat = (122 5—)0.5*ziovL )**.333333 
phi2 = phil 
phis = ( l. = 14.*zovL )**.25 
Pilop= cl. + .75*(—zovL) **.6667)**1.5 
epsilonO = ustar3*phiEp*zinv 
GOTO 440 
430 CONTINUE 


kkk kkk kkkkK 


c¢ get Km, Ri, and other parameters for neutral to stable BL 
kKkkkkkkk kk 
h = zi 
IF ( zi .lt. -98. ) h = sqrt(k*ustar*Linv/f) 
a IF(i.lt.iw)WRITE(6,*)'h=zi,or Vv (ku*/(fL))',h,zi,f 
hinv = 1./h 
ZOVH = z*hinv 
ZSQOVHL = zovL*zovH 
ZSQOVH2 = zovH*zovH 
alpha = 2.0 - 10.0*Linv 
beta = 3.0 - 20.0*Linv 
d IF(i.lt.iw)WRITE(6,*)‘'a,8,z/H',alpha, beta, zovH 
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tps = 1,.-—. 2ovil 
tp8alpha = tp8**alpha 
a IF(i.1t.iw)WRITE(6,*)'(1-z/H),(1-z/H)**a", tps, epsalpna 
tp8al2 = sqrt(tp8alpha) 
tp8beta = tp8&**beta 


d IF(i.1lt.iw)WRITE(6,*)'(1-z/H)**a/2,(1-z/H)**B',tp8al2,tp8beta 
phil = 2.2*tp8al2 
phi2 = phil 
phi3 = 1.6*tp8al2 

e local friction velocity 


ul = ustar*sqrt (tp8alpha) 
IF(i.l1t.iw)WRITE(6,*) ‘ul=u*(1-z/H)**(a/2)',ul,ustar, zovH 
local Obukhov length 
Ll = L*tp8**(1.5*alpha - beta) 
IF(i.lt.iw)WRITE(6,*) 'L1=L(1-z/H) **(3a/2-8)',L1,L,zovH, alpha, 
& beta 
total shear 
adUtotdz = 4.7*ul/(k*Ll) 
dUtotdz = 2.5*ustar*zinv*(1.+4.7*zovL) **(0.5*alpha) *tp8al2 
IF(i.1lt.iw)WRITE(6,*) 'dU/dz=4.7ul/(kL1)',dUtotdz,ul,k,Ll 
temperature and buoyancy fluxes at height, z 
wpTpZ = wpTp0*(1l. - zovzi)**beta 
Brunt Vaisala frequency at z 
tp9 = (1. + 3.7*ZovE) <ziny 
IF(i.lt.iw)WRITE(6,*) 'tp9=(1+3.7z/L)/z',tp9,zovL, zinv 
BVF = 4.3*ustar*Tstar**2*tp9*tp8beta*tp8alpha 
IF(i.lt.iw)WRITE(6,*) 'BVF=4. 3u*O*2(1+3.72/L) (1-z/H) **(a+B)/z', 
& BVF,ustar,Tstar, ZovH, ZOvVL,Z 
TKE dissipation rate at z 
phiEp = 3.6*tp9*tp8beta*zinv 
epsilon = phiEp*ustar3 
c Richardson number at z 
tploO = 1./(1.. + S.*zovL.) 
Ri = zovL*tp10 


Q Qo qgaag aa Q a 


aaa O. 


e momentum/heat diffusivities at z 
Km = k*ustar*z*tp8*tplo 
a IF(i.lt.iw)WRITE(6,*) 'Km=ku*z(1-z/H) /(1+5z/L)',Km,k,ustar,z, 
a & ZOVH, ZOVL 
KkKAKKKKKKKK KKK KKK KK KKK KR KKK 
440 CONTINUE 


c Get ou,v,w, tke, total horizontal velocity, & temperature 
SigmaU = ustar*phil 
SigmaU2= sigmaU*sigmaU 
SigmaV = 0.894*sigmaU 
SigmaV2= 0.5*sigmaU2 
SigmaW = ustar*phi3 
SigmaW2 = sigmaW**2 
tke = 0.5*(sSigmaU2 + sigmaV2 + sigmaW2) 
d IF(i.1t.iw)WRITE(6,*) 'tke=oui2/2',tke,sigmaU2,sigmaV2,sigmaW2 
Cc IF ( icase .eq. 2 ) tke = tkeusl 
Oldtheta = Theta 
Theta = Theta2 + 2.5*Tstar*( lnzovzO - psih ) 
IF ( i. eq. 1) Oldtheta = Theta 
Deltheta = Theta —- Oldtheta 


da IF(i.lt.iw)WRITE (6,*)'50,0,01d0', Deltheta, Theta, Oldtheta 
460 CONTINUE 
Vtotal = 2.5*ustar*(lnzovzO - psim) 
d IF(i.lt.iw)WRITE(6,*) 'V=2.5u*(Ln(z/z0)-psim))',Vtotal,ustar, 
a & lnzovz0O,psim 
u = Vtotal 
a IF(i.1lt.iw)WRITE(6,*) '9=02+2.50*(Ln(z/z0)-psih)',theta,theta2, 
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re | & tstar, lnzovz0,psih 
KKK KKK KK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KEK KKK KKK KEKE KKK KKEKKEKKEKE 
c debug formats 
d1000 FORMAT (10x, 2911.3) 
@1010 FORMAT (10x, 3f11.3) 
d1020 FORMAT (10x, 4f11.3) 
a1030 FORMAT (1x, i5 ) 
d1040 FORMAT (10x, 3g11.3, i5) 
ad1050 FORMAT (10x, 4f11.3, i5) 
eueos FORMAT (10x, Sf11.3, i5) 
@a1060 CONTINUE 
RETURN 
END 


KEKE KKEKKKKEKK KEK KKK KE KKK KEKE KKK KKK KKK KKEKKEKKEKEKKEKKKEKKKEKEKEKKEKKEKKKEKEKKEKKKEKKKKK 


SUBROUTINE NORNG(r,p) 


G 
c This subroutine generates a sequence of numbers normally and randomly 
c distributed over the interval -3 to 3 from uniformly distributed random 
c numbers, by the method of linear approximation to the inverse of the 
c¢ accumulative normal distribution function 
= 
DOUBLE PRECISION r,p, y(6), x(6), 3(5) 
BATA y/O7d0, 0.0228d0, ©.0665d0, 0.1357d0, 0.2743d0, 0.5d0 /, 
& x/ -3.01d0, -270d0, -1.5d07">-1:0d0, -0.6d0, 0.0da0 /, 
& s/ 43.8596d0, 11.3636d0, 7.25689d0, 2.891352d0, 2.65887d0 / 
CALL STRNUM(r) 
p=r 
i=l 
ie (spe. gt 0.500.) p= 1.0d0"=— r2 EE ( eett.s ytitl) ) GOTO:s 
i=it+l 
GOTO 2 
8 P= ((p- yi) )*s(i) + x(i) ) 
tier .ge. O.5d0 ) p = =p 
cd WRITE (6,*)'random numbers, r & p',r,p 
RETURN 
END 


REKKEKEKEKKEKKKEKKKEKEKKEKEKEKKRKKKKKKKKKKKKEKKKKKEKKKKKEKKKKEKEKRKKKEKKKEKKKKKKKKKKKKKKEKKK 


SUBROUTINE STRNUM(r) 


Cc * 
c This subroutine generates a sequence of numbers which are randomly and 
c uniformly distributed over the unit interval. 
c 

DOUBLE PRECISION pl, r 

bb = 1.d0 

pl = r*317.d0 
c WRITE (6, *)*pl’,pl 

r = abs( pl - ( INT(pl/bb)*bb) ) 

RETURN 

END 


KEKE KKKEKKKEKKKEKEKKKEKEK KKK KEK KEK KKK KKKKEKEKKEKEKEKEKKKEKKEKKKKKKKKKKKKK KK KKK KKKKKKK 


FUNCTION MAX (a, b) 
max = b 
IF ( a .gt. b) max =a 
END 
KKK KKKKKEKEK KEK KKK KK KKK KKK KEKE KEK KKEKKKKKKEKK KKK KKK KKK KK KKK KKK KKK KK KKK KKEKK 
FUNCTION MIN (a, b) 
min = b 
IF ( a .lt. b) max =a 
END 


RKKKEKKRKKKRK KKK KK KEKE KKK K KKK KEKE KKK KEKE KKEKKEKEKEEKEKEKKKKEKEKREKKKKKEKKKEKEKEKKEKKKKKRKRKKRKEK 
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